### Replication Package for "Why is Intermediating Houses so Difficult? Evidence from iBuyers"
### Buchak, Matvos, Piskorski, and Seru
###
###
### buchak@stanford.edu


library(data.table)
library(ggplot2)
library(nleqslv)
library(stats)
library(doParallel)
library(stargazer)


cl = makeCluster(5)
registerDoParallel(cl)

# Toggle to reestimate model
REESTIMATE_MODEL = T

# Toggle to create counterfactuals
RECREATE_COUNTERFACTUALS = T


run_all <- function() {

  if(REESTIMATE_MODEL) {
    calibrateModel()
  }
  
  createOutput(RECREATE_COUNTERFACTUALS)
  
}







calibrateModel <- function() {
  
  
  gmmObjectiveALL <- function(x,return_solved_model = F) {
    
    # Deal parameters
    new.pars <- rbind(
      data.table(variable = 'u_ib_delta'       ,value = x[1]),
      data.table(variable = 'lambda_ib'        ,value = x[2]),
      data.table(variable = 'ch'               ,value = x[3]),
      data.table(variable = 'nu'               ,value = x[4]),
      data.table(variable = 'sig_i'            ,value = x[5]),
      data.table(variable = 'signal_noise'     ,value = x[6]),
      data.table(variable = 'u_pat_delta'   ,value = x[7]),
      data.table(variable = 'u_imp_delta'   ,value = x[8]),
      data.table(variable = 'u_delta_bad'   ,value = x[9]),
      data.table(variable = 'p_switch_types',value = x[10]),
      data.table(variable = 'lambda'        ,value = x[11]),
      data.table(variable = 'sig_m',         value = x[12]),
      data.table(variable = 'ibu',           value = x[13]))
    
    
    print('---------------------------- PARAMETERS ----------------------------')
    print(new.pars)
    print('--------------------------------------------------------------------')
    print('')
    
    
    all.pars <- new.pars
    pars <- setPars(all.pars)
    
    # Solve model
    solved.model      <- solveModelRevision(pars = pars,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = ibuyer.guess,pre_ibuyer_only = F)
    share.derivatives <- getIBShareDerivatives(all.pars,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = ibuyer.guess)
    
    # Get moments
    model.moments <- calculateiBuyerMoments(solved.model,share.derivatives)
    
    # Merge moments
    merged.moments <- merge(ibuyer.moments,model.moments,by='name',suffixes=c('.empirical','.model'))
    
    
    # Calculate GMM objective
    merged.moments[,dist := value.empirical - value.model]
    merged.moments[,d2   := dist^2]
    merged.moments[,cont := d2 * weight]
    merged.moments[is.na(cont),cont := 99999999999999999]
    
    # Add line for positive profitability
    profitability   = sum(solved.model$ibuyer$stats$ib_prod_share$share * solved.model$ibuyer$stats$ib_prod_share$profitability)/solved.model$no.ibuyer$stats$price_hh_mean
    
    
    merged.moments <- rbind(merged.moments,data.table(name = 'profitability',
                                                      value.empirical = 0,
                                                      value.model = profitability,
                                                      weight=1,
                                                      dist = (profitability < 0) * abs(profitability),
                                                      d2 = ((profitability < 0) * abs(profitability))^2,
                                                      cont = 1000000 * ((profitability < 0) * abs(profitability))^2))
    
    
    gmm.objective = sum(merged.moments$cont)
    
    merged.moments <- merged.moments[order(-cont)]
    
    # Update initial guess if we have a better one
    if(gmm.objective < .95 * best_gmm_objective) {
      best_gmm_objective <<- gmm.objective
      ibuyer.guess <<- solved.model$ibuyer$solution
    }
    
    # Print some status updates
    print('----------------------------- MOMENTS ------------------------------')
    print(merged.moments[,c('name','value.empirical','value.model','cont')])
    print('--------------------------------------------------------------------')
    print('')
    print('---------------------------- OBJECTIVE -----------------------------') 
    print(gmm.objective)
    print('--------------------------------------------------------------------')
    print('')
    print('')
    print('')
    
    aa <- solved.model$ibuyer$stats$ib_prod_share
    aa <- aa[order(-share)]
    print(aa[1:10])
    
    if(!return_solved_model) {
      return(gmm.objective)
    } else {
      return(list(moments = merged.moments,objective = gmm.objective,params = pars,par.dt = all.pars,solved.model = solved.model,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = ibuyer.guess))
    }
    
  }
  
  
  X = 1
  ibuyer.moments <- getEmpiricalMoments_Post_iBuyer()
  pre.ibuyer.moments <- getEmpiricalMoments_Pre_iBuyer()
  
  
  no.ibuyer.guess <<- c(317.6541000, 311.0227000, 234.4917000, 227.2361000, 
                        307.2213000, 299.1472000, 224.2220000, 215.5249000,   
                        0.0000000,   0.4659618,   0.4646781)
  
  
  ibuyer.guess <- c(3.048009e+02, 2.975219e+02, 2.545778e+02, 2.469231e+02,
                    2.861982e+02, 2.376091e+02, 2.374365e+02, 2.374345e+02,
                    3.056280e+02, 2.908273e+02, 2.552824e+02, 2.400673e+02,
                    2.936322e+02, 2.847970e+02, 2.435073e+02, 2.342972e+02,
                    4.662764e-01, 4.653042e-01,
                    5.270738e-04, 9.386437e-05, 1.069467e-03, 1.925519e-04)
  
  
  best_gmm_objective <<- 99999999
  
  guess   =      c( 0.26681042 , 0.97930652,146.77819994,  1.66325260 ,  6.35775674 ,  0.03351957, 36.20798584  ,40.14513868   ,6.20747437        ,0.24345591  ,724.21998090  ,  9.77496571,-14.76126797)
  scale   =      1e-2*c(1   ,.1  ,5       ,.1    ,.1 , .15,1  ,1   ,.1       , .1    ,10   ,.1,2)
  bound_l =      c(-10 ,.975 ,.01        ,.220  , .75  , .005,10  ,10   ,.75        ,.015   ,10    ,1,-50) 
  bound_u =      c(150 ,1.05 ,350      ,50     ,15  , .135,50  ,50  ,25       ,0.35    ,750  ,25,25)

  
  rr1 <- optim(par = guess, fn = gmmObjectiveALL,control = list(maxit = 1000,parscale = scale,lmm = 4),method = 'L-BFGS-B',lower = bound_l,upper = bound_u)

  
  calibrated.model <- gmmObjectiveALL(rr1$par,return_solved_model = T)

  
  save(calibrated.model,file='estimates_and_counterfactuals/calibration_production.rsave')

}



createOutput <- function(RECREATE_COUNTERFACTUALS = F) {
  
  # Load calibrated parameters
  load('estimates_and_counterfactuals/calibration_production.rsave') 
  
  # 0. Table 5: Moments, parameters
  
  # Table 5 Panel A (Moments)
  moments <- calibrated.model$moments[,c('name','value.empirical','value.model')]
  name.dt <- data.table(name   = c('share.unprofitable','buy.discount','p.reduce.prices','diB.dErr','ibuyer.spread','price_sd_xs','diB.dLiq','imp.price.delta.pct','iBuyer.share','tom.hh','price_sd_ts','tom.ib.delta','imp.list.delta.days','price_listed'),
                        Moment = c('Share unprofitable (%)','iBuyer buy discount (%)','P(iBuyer cut prices) (%)','d(iB share)/d(Pricing error)','iBuyer gross spread (%)','XS price variation (%)','d(iB share)/d(P(sells in 90 days))','Impatient price delta (%)','iBuyer market share (%)','HH time on market (days)','TS price variation (%)','iBuyer time on market (days)','Impatient time on market delta (days)','List price (k)'))
  moments <- merge(moments,name.dt,by='name')[,c('Moment','value.empirical','value.model')]
  names(moments) <- c('Moment','Data','Model')
  moments <- moments[c(9,13,11,10,7,6,4,14,5,1,12,3,2,8)]
  
  # Moment is actually iBuyer DELTA time on market but we display the total time on market.
  moments[Moment == 'iBuyer time on market (days)']$Data <-  moments[Moment == 'iBuyer time on market (days)']$Data + moments[Moment == 'HH time on market (days)']$Data
  moments[Moment == 'iBuyer time on market (days)']$Model <-  moments[Moment == 'iBuyer time on market (days)']$Model + moments[Moment == 'HH time on market (days)']$Model
  
  moments[c(3,4,5,7,9,10,11,14)]$Data <- moments[c(3,4,5,7,9,10,11,14)]$Data*100 # Display as percentages
  moments[c(3,4,5,7,9,10,11,14)]$Model <- moments[c(3,4,5,7,9,10,11,14)]$Model*100  # Display as percentages
  moments[Moment == 'iBuyer buy discount (%)']$Model <- moments[Moment == 'iBuyer buy discount (%)']$Model * -1  # Quote as discount (rather than negative premium)
  moments[Moment == 'iBuyer buy discount (%)']$Model <- moments[Moment == 'iBuyer buy discount (%)']$Model * -1  # Quote as discount (rather than negative premium)
  
  
  stargazer(moments,summary=F,out = '../out/tables/5A.html',type = 'html')
  
  # Table 5 Panel B Is manually created (see sources in table).
  
  # Table 5 Panel C (Calibrated parameters)
  par.dt <- calibrated.model$par.dt
  par.dt$variable <- c('m_l','lambda_ib','mh_minus_ml','eta','sigma_i','xi','delta_u_patient','delta_u_impatient','u_bad','p_switch','lambda','sigma_m','delta_ib')
  par.dt <- par.dt[c(7,8,13,1,4,3,11,2,6,9,12,5,10)]
  par.dt[3]$value <- par.dt[3]$value*-1 # Quoted as disutility rather than negative utility
  stargazer(par.dt,summary=F,out = '../out/tables/5C.html',type = 'html')
  
  
  
  # Table 6 Panels A, B, and C: Compare ibuyer vs. non-ibuyer equilibrium
  no.ib <- calibrated.model$solved.model$no.ibuyer
  ib    <- calibrated.model$solved.model$ibuyer
  comparison = rbind(data.table(variable = 'avg_price',no.ib = no.ib$stats$price_hh_mean,ib = ib$stats$ib_share * ib$stats$price_ib_buy_mean + (1-ib$stats$ib_share)*ib$stats$price_hh_mean),
                     data.table(variable = 'hh_price',no.ib = no.ib$stats$price_hh_mean,ib = ib$stats$price_hh_mean),
                     data.table(variable = 'tom_avg' ,no.ib = no.ib$stats$tom_hh_mean,ib = ib$stats$ib_share * 0 + (1-ib$stats$ib_share)*ib$stats$tom_hh_mean),
                     data.table(variable = 'tom_hh'  ,no.ib = no.ib$stats$tom_hh_mean,ib = ib$stats$tom_hh_mean),
                     data.table(variable = 'v_hh'    ,no.ib = no.ib$stats$vhh,ib = ib$stats$vhh),
                     data.table(variable = 'v_h_g'   ,no.ib = no.ib$stats$vh_good,ib=ib$stats$vh_good),
                     data.table(variable = 'v_h_b'   ,no.ib = no.ib$stats$vh_bad,ib=ib$stats$vh_bad),
                     data.table(variable = 'v_s_pg'   ,no.ib = no.ib$stats$vs_good_pat,ib=ib$stats$vs_good_pat),
                     data.table(variable = 'v_s_pb'   ,no.ib = no.ib$stats$vs_bad_pat,ib=ib$stats$vs_bad_pat),
                     data.table(variable = 'v_s_ig'   ,no.ib = no.ib$stats$vs_good_imp,ib=ib$stats$vs_good_imp),
                     data.table(variable = 'v_s_ib'   ,no.ib = no.ib$stats$vs_bad_pat,ib=ib$stats$vs_bad_pat))
  comparison[,pct.change := ib / no.ib-1]                   
  
  p6a <- comparison[variable %in% c('avg_price','hh_price','tom_avg','tom_hh'),c('variable','pct.change')]
  p6a$pct.change <- p6a$pct.change*100 # quote in percent
  p6a$variable <- c('Price (all transactions)','Price (h2h)','Time on market (all transactions)','Time on market (h2h)')
  stargazer(p6a,summary=F,out='../out/tables/6A.html',type = 'html')
  
  p6b <- comparison[variable %in% c('v_hh','v_h_g','v_h_b','v_s_pg','v_s_pb','v_s_ig','v_s_ib'),c('variable','pct.change')]
  p6b$pct.change <- p6b$pct.change * 100
  p6b$variable <- c('Unconditional average','Matched (good house)','Matched (bad house)','Patient seller (good house)','Patient seller (bad house)','Impatient seller (good house)','Impatient seller (bad house)')
  p6b <- rbind(p6b,data.table(variable = 'Buyer',pct.change = 0))
  stargazer(p6b,summary=F,out='../out/tables/6B.html',type = 'html')
  
  
  # Table 6 Panel C: Shares and Time on Market by Seller Type
  comparison = rbind(data.table(variable = 'iBuyer share of patient',no.ib = 0,ib = sum(ib$stats$ib_prod_share[patient == 'p']$share_total/.50)), # .5 of total are patient
                     data.table(variable = 'iBuyer share of impatient',no.ib = 0,ib = sum(ib$stats$ib_prod_share[patient == 'i']$share_total/.50)), # .5 of total are impatient
                     data.table(variable = 'Impatient share of iBuyer',no.ib = 0,ib = sum(ib$stats$ib_prod_share[patient == 'i']$share)),
                     data.table(variable = 'Time on market patient',no.ib = no.ib$stats$tom_hh_patient_mean,ib = ib$stats$tom_hh_patient_mean * (1-sum(ib$stats$ib_prod_share[patient == 'p']$share_total/.50))),
                     data.table(variable = 'Time on market impatient',no.ib = no.ib$stats$tom_hh_impatient_mean,ib = ib$stats$tom_hh_impatient_mean * (1-sum(ib$stats$ib_prod_share[patient == 'i']$share_total/.50))))
  comparison[,Change := ib / no.ib - 1]
  comparison[is.infinite(Change),Change := 0]
  comparison[variable %in% c('iBuyer share of patient','iBuyer share of impatient','Impatient share of iBuyer')]$ib <- 100 * comparison[variable %in% c('iBuyer share of patient','iBuyer share of impatient','Impatient share of iBuyer')]$ib
  comparison[variable %in% c('Time on market patient','Time on market impatient')]$Change <- 100 * comparison[variable %in% c('Time on market patient','Time on market impatient')]$Change
  names(comparison) <- c('Outcome','No-iBuyer equilibrium','iBuyer Equilibrium','% Change')
  
  stargazer(comparison,summary=F,out='../out/tables/6C.html',type = 'html')
  
  
  
  if(RECREATE_COUNTERFACTUALS) {
  
    N = 1000
    
    
    # 1.a Run comparative statics.
    liquidity.cs    <- runComparativeStatic(par = 'lambda',lower = .05, upper = 50, length = N,calibrated.model)
    signal.cs       <- runComparativeStatic(par = 'signal_noise',lower = .1, upper = 5, length = N*.5,calibrated.model) 
    delay.cs        <- runComparativeStatic(par = 'delay',lower = .25, upper = 3, length = N,calibrated.model)
    rent.cs        <- runComparativeStatic(par = c('u_ib_delta','ch'),lower = 0, upper = 1.1, length = N,calibrated.model)
    dispersion.cs        <- runComparativeStatic(par = 'dispersion',lower = -3, upper = 7, length = N,calibrated.model)
    preference.cs  <- runComparativeStatic(par = 'ibu',lower = 0,upper = 2,length=N,calibrated.model)
    heterogeneity.cs  <- runComparativeStatic(par = 'sig_m',lower = .25,upper = 1.75,length=N,calibrated.model)
    match.cs      <- runComparativeStatic(par = 'lambda_ib',lower = .5,upper = 2,length=N,calibrated.model)
    
    # Do one where the iBuyer has no noise but is slow
    no.noise.model <- calibrated.model
    no.noise.model$par.dt[variable == 'signal_noise']$value <- 0
    no.noise.delay.cs        <- runComparativeStatic(par = 'delay',lower = .25, upper = 3, length = N,no.noise.model)
    
    
    save(list=c('no.noise.delay.cs','liquidity.cs','signal.cs','delay.cs','rent.cs','dispersion.cs','heterogeneity.cs','match.cs','preference.cs'),file = 'estimates_and_counterfactuals/saved_counterfactuals_baseline.rsave')
    
    # 1.b. Run comparative statics as if u_ib_delta is 0
    # Model with no anti-ib preference Cor some (appendix) counterfactuals
    no.disutility.model = calibrated.model
    no.disutility.model$par.dt[variable == 'ibu']$value <- 0
    
    liquidity.cs    <- runComparativeStatic(par = 'lambda',lower = .025, upper = 100, length = N,no.disutility.model)
    signal.cs       <- runComparativeStatic(par = 'signal_noise',lower = .05, upper = 15, length = N,no.disutility.model)
    delay.cs        <- runComparativeStatic(par = 'delay',lower = .25, upper = 10, length = N,no.disutility.model)
    rent.cs        <- runComparativeStatic(par = c('u_ib_delta','ch'),lower = 0, upper = 1.1, length = N,no.disutility.model)
    dispersion.cs        <- runComparativeStatic(par = 'dispersion',lower = -5, upper = 10, length = N,no.disutility.model)
    preference.cs  <- runComparativeStatic(par = 'ibu',lower = 0,upper = 3,length=N,no.disutility.model)
    heterogeneity.cs  <- runComparativeStatic(par = 'sig_m',lower = .05,upper = 2.75,length=N,no.disutility.model)
    match.cs      <- runComparativeStatic(par = 'lambda_ib',lower = .5,upper = 2,length=N,no.disutility.model)
    
    save(list=c('liquidity.cs','signal.cs','delay.cs','rent.cs','dispersion.cs','heterogeneity.cs','match.cs'),file = 'estimates_and_counterfactuals/saved_counterfactuals_no_hedonic_disutility.rsave')
    
    
    
  }
  
  
  # Load counterfactuals and create charts
  load('estimates_and_counterfactuals/saved_counterfactuals_baseline.rsave')
  
  
  ###### Figure 2: One-off decompositions (side-by-side) #####
  # 
  # 1. Baseline
  # 2. Slow (delay ~= 45/360)
  # 3. Rent (rent params = 0)
  # 4. Signal (signal = .125)
  speed.target = 30/360
  rent.target = rent.cs$param.value[1] * (1-.5) # https://www.zillow.com/learn/investing-101-estimating-rental-property-expenses/ 0.65 means the costs are 35% of gross income
  signal.target = .5
  
  liquidity.cs[,d :=  0]
  delay.cs[,d :=  (param.value - speed.target)^2]
  rent.cs[,d :=   (param.value - rent.target)^2]
  signal.cs[,d := (param.value - signal.target)^2]
  no.noise.delay.cs[,d :=  (param.value - speed.target)^2]
  
  baseline.cf <- liquidity.cs[scale == 1][1]
  speed.cf    <- delay.cs[d == min(delay.cs$d)][1]
  rent.cf     <- rent.cs[d == min(rent.cs$d)][1]
  signal.cf   <- signal.cs[d == min(signal.cs$d)][1]
  no.nose.speed.cf <- no.noise.delay.cs[d == min(delay.cs$d)][1]
  
  baseline.cf[,Scenario := 'Baseline']
  speed.cf[,Scenario := 'Slow\niBuyer']
  rent.cf[,Scenario := 'Renting\niBuyer']
  signal.cf[,Scenario := 'Inaccurate\niBuyer']
  no.nose.speed.cf[,Scenario := 'No noise/slow\niBuyer']
  
  hh = 4
  ww = 4
  
  combined <- rbind(baseline.cf,speed.cf,rent.cf,signal.cf)
  combined[,Scenario := factor(Scenario,c('Baseline','Inaccurate\niBuyer','Slow\niBuyer','Renting\niBuyer'))]
  
  
  combined2 <- rbind(baseline.cf,no.nose.speed.cf)
  combined2[,Scenario := factor(Scenario,c('Baseline','Inaccurate\niBuyer','No noise/slow\niBuyer','Slow\niBuyer','Renting\niBuyer'))]
  
  ggplot(combined) + geom_bar(aes(x=Scenario,y=100*iBuyer_share),color='black',stat = 'identity',position = 'dodge') + ylab('%') + xlab(NULL) + theme_bw() + 
    theme(legend.position = 'bottom',legend.title = element_blank())
  ggsave('../out/figures/2A.png',height=hh,width=ww,units='in',dpi = 300)
  
  ggplot(combined) + geom_bar(aes(x=Scenario,y=100*profitability),color='black',stat = 'identity',position = 'dodge') + ylab('%') + xlab(NULL) + theme_bw() + 
    theme(legend.position = 'bottom',legend.title = element_blank())
  ggsave('../out/figures/2B.png',height=hh,width=ww,units='in',dpi = 300)
  
  ggplot(combined) + geom_bar(aes(x=Scenario,y=100*share_imp),color='black',stat = 'identity',position = 'dodge') + ylab('%') + xlab(NULL) + theme_bw() + 
    theme(legend.position = 'bottom',legend.title = element_blank())
  ggsave('../out/figures/2C.png',height=hh,width=ww,units='in',dpi = 300)
  
  ggplot(combined) + geom_bar(aes(x=Scenario,y=100*share_b),color='black',stat = 'identity',position = 'dodge') + ylab('%') + xlab(NULL) + theme_bw() + 
    theme(legend.position = 'bottom',legend.title = element_blank())
  ggsave('../out/figures/2D.png',height=hh,width=ww,units='in',dpi = 300)
  
  
  
  
  ##### Continuous versions of these 
  
  # Clean up discontinuities and unsolved equilibria
  liquidity.cs[profitability < 0 | dist > 1e-8,iBuyer_share := NA]
  liquidity.cs[profitability < 0 | dist > 1e-8,share_b := NA]
  liquidity.cs[profitability < 0 | dist > 1e-8,share_imp := NA]
  liquidity.cs[profitability < 0 | dist > 1e-8,profitability := NA]
  
  signal.cs[profitability < 0 | dist > 1e-8,iBuyer_share := NA]
  signal.cs[profitability < 0 | dist > 1e-8,share_b := NA]
  signal.cs[profitability < 0 | dist > 1e-8,share_imp := NA]
  signal.cs[profitability < 0 | dist > 1e-8,profitability := NA]
  
  signal.cs <- signal.cs[order(param.value)]
  signal.cs[,dShare := iBuyer_share - shift(iBuyer_share)]
  signal.cs[dShare < -.01,iBuyer_share := NaN]
  signal.cs[dShare < -.01,profitability := NaN]
  signal.cs[dShare < -.01,share_imp := NaN]
  signal.cs[dShare < -.01,share_b := NaN]
  signal.cs[dShare < -.01,house_price_ib := NaN]
  
  dispersion.cs[profitability < 0 | dist > 1e-8,iBuyer_share := 0]
  dispersion.cs[profitability < 0 | dist > 1e-8,share_b := NA]
  dispersion.cs[profitability < 0 | dist > 1e-8,share_imp := NA]
  dispersion.cs[profitability < 0 | dist > 1e-8,profitability := 0]
  
  # For signal, if there's a big drop, we NaN it out.
  signal.cs <- signal.cs[order(param.value)]
  signal.cs[,dShare := iBuyer_share - data.table::shift(iBuyer_share)]
  signal.cs[abs(dShare) > 0.05,iBuyer_share := NaN]
  signal.cs[abs(dShare) > 0.05,profitability := NaN]
  signal.cs[abs(dShare) > 0.05,share_imp := NaN]
  signal.cs[abs(dShare) > 0.05,share_b := NaN]
  
  signal.cs[,dHP := house_price_ib/shift(house_price_ib) - 1]
  signal.cs[,dW := welfare_nib/shift(welfare_nib) - 1]
  
  signal.cs[abs(dHP) > 0.05,house_price_ib := NaN]
  signal.cs[abs(dW) > 0.0008,welfare_nib := NaN]
  
  
  signal.cs <- signal.cs[param.value <= 0.75 ]
  
  # Signal
  ggplot(signal.cs) + geom_line(aes(x=param.value,y=100*iBuyer_share)) + theme_bw() + xlab('iBuyer signal noise') + ylab('%') + geom_vline(xintercept = calibrated.model$par.dt[variable == 'signal_noise']$value,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/3A.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(signal.cs) + geom_line(aes(x=param.value,y=100*profitability)) + theme_bw() + xlab('iBuyer signal noise') + ylab('%') + geom_vline(xintercept = calibrated.model$par.dt[variable == 'signal_noise']$value,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/3B.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(signal.cs) + geom_line(aes(x=param.value,y=100*share_imp)) + theme_bw() + xlab('iBuyer signal noise') + ylab('%') + geom_vline(xintercept = calibrated.model$par.dt[variable == 'signal_noise']$value,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/3C.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(signal.cs) + geom_line(aes(x=param.value,y=100*share_b)) + theme_bw() + xlab('iBuyer signal noise') + ylab('%') + geom_vline(xintercept = calibrated.model$par.dt[variable == 'signal_noise']$value,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/3D.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  # Remove big spikes in house price due to numerical issues
  liquidity.cs <- liquidity.cs[order(tom)]
  liquidity.cs[,dHP := house_price_ib/house_price_nib - data.table::shift(house_price_ib/house_price_nib) ]
  liquidity.cs[,dShare := iBuyer_share - data.table::shift(iBuyer_share)]
  liquidity.cs[abs(dShare)>0.0009,iBuyer_share := NaN]
  liquidity.cs[abs(dShare)>0.0009,profitability := NaN]
  liquidity.cs[abs(dShare)>0.0009,share_imp := NaN]
  liquidity.cs[abs(dShare)>0.0009,share_b := NaN]
  liquidity.cs[abs(dShare)>0.0009,house_price_ib := NaN]
  
  
  #liquidity.cs <- liquidity.cs[abs(dHP)<0.11 & (house_price_ib/house_price_nib-1) < .3]
  #liquidity.cs[abs(dHP)>0.001,house_price_ib := NaN]
  
  # Liquidity
  ggplot(liquidity.cs) + geom_line(aes(x=tom,y=100*iBuyer_share)) + theme_bw() + xlab('Time on market') + ylab('%') + geom_vline(xintercept = baseline.cf$tom,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/4A.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(liquidity.cs) + geom_line(aes(x=tom,y=100*profitability)) + theme_bw() + xlab('Time on market') + ylab('%') + geom_vline(xintercept = baseline.cf$tom,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/4B.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(liquidity.cs) + geom_line(aes(x=tom,y=100*share_imp)) + theme_bw() + xlab('Time on market') + ylab('%') + geom_vline(xintercept = baseline.cf$tom,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/4C.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(liquidity.cs) + geom_line(aes(x=tom,y=100*share_b)) + theme_bw() + xlab('Time on market') + ylab('%') + geom_vline(xintercept = baseline.cf$tom,linetype = 'dashed')  + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/4D.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  
  heterogeneity.cs <- heterogeneity.cs[house_price_ib/house_price_nib-1 > 0 & scale > 0  & (house_price_ib/house_price_nib-1) < .25]
  heterogeneity.cs <- heterogeneity.cs[order(scale)]
  ggplot(heterogeneity.cs) + geom_line(aes(x=scale,y=100*iBuyer_share)) + theme_bw() + xlab('Dispersion of private values') + ylab('%') + geom_vline(xintercept = 1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/5A.png',height=hh,width=ww,units = 'in',dpi = 300)
  
    ggplot(heterogeneity.cs) + geom_line(aes(x=scale,y=100*profitability)) + theme_bw() + xlab('Dispersion of private values') + ylab('%') + geom_vline(xintercept = 1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/5B.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(heterogeneity.cs) + geom_line(aes(x=scale,y=100*share_imp)) + theme_bw() + xlab('Dispersion of private values') + ylab('%') + geom_vline(xintercept = 1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/5C.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(heterogeneity.cs) + geom_line(aes(x=scale,y=100*share_b)) + theme_bw() + xlab('Dispersion of private values') + ylab('%') + geom_vline(xintercept = 1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/5D.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  
  # Appendix Figure A11: Changing Hedonic Preference
  preference.cs.plot <- preference.cs[scale <= 1.1]
  ggplot(preference.cs.plot) + geom_line(aes(x=-1*scale,y=100*iBuyer_share)) + theme_bw() + xlab('Hedonic disutility from iBuyer sale') + ylab('%') + geom_vline(xintercept = -1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/A11A.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(preference.cs.plot) + geom_line(aes(x=-1*scale,y=100*profitability)) + theme_bw() + xlab('Hedonic disutility from iBuyer sale') + ylab('%') + geom_vline(xintercept = -1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/A11B.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(preference.cs.plot) + geom_line(aes(x=-1*scale,y=100*share_imp)) + theme_bw() + xlab('Hedonic disutility from iBuyer sale') + ylab('%') + geom_vline(xintercept = -1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/A11C.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  ggplot(preference.cs.plot) + geom_line(aes(x=-1*scale,y=100*share_b)) + theme_bw() + xlab('Hedonic disutility from iBuyer sale') + ylab('%') + geom_vline(xintercept = -1,linetype = 'dashed') + geom_hline(yintercept = 0,linetype = 'dashed')
  ggsave('../out/figures/A11D.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  # Appendix Figure A12: Some comparative statics with hedonic disutility set to 0
  load('estimates_and_counterfactuals/saved_counterfactuals_no_hedonic_disutility.rsave')
  # Clean up discontinuities and unsolved equilibria
  liquidity.cs[profitability < 0 | dist > 1e-8,iBuyer_share := NA]
  signal.cs[profitability < 0 | dist > 1e-8,iBuyer_share := NA]
  dispersion.cs[profitability < 0 | dist > 1e-8,iBuyer_share := 0]
  
  # For signal, if there's a big drop, we NaN it out.
  signal.cs <- signal.cs[order(param.value)]
  signal.cs[,dShare := iBuyer_share - data.table::shift(iBuyer_share)]
  signal.cs[abs(dShare) > 0.05,iBuyer_share := NaN]
  
  signal.cs <- signal.cs[param.value <= 0.75 & param.value > 0.01]
  
  # Signal
  ggplot(signal.cs) + geom_line(aes(x=param.value,y=100*iBuyer_share)) + theme_bw() + xlab('iBuyer signal noise') + ylab('%') + geom_vline(xintercept = calibrated.model$par.dt[variable == 'signal_noise']$value,linetype = 'dashed')  
  ggsave('../out/figures/A12A.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  # Remove big spikes in house price due to numerical issues
  liquidity.cs <- liquidity.cs[order(tom)]
  liquidity.cs[,dShare := iBuyer_share - data.table::shift(iBuyer_share)]
  liquidity.cs <- liquidity.cs[abs(dShare) < 0.006 | tom > 100]
  
    # Liquidity
  ggplot(liquidity.cs) + geom_line(aes(x=tom,y=100*iBuyer_share)) + theme_bw() + xlab('Time on market') + ylab('%') + geom_vline(xintercept = baseline.cf$tom,linetype = 'dashed') 
  ggsave('../out/figures/A12B.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  # House heterogeneity
  heterogeneity.cs <- heterogeneity.cs[house_price_ib/house_price_nib-1 > 0 & scale > 0  & (house_price_ib/house_price_nib-1) < .25]
  heterogeneity.cs <- heterogeneity.cs[order(scale)]
  ggplot(heterogeneity.cs) + geom_line(aes(x=scale,y=100*iBuyer_share)) + theme_bw() + xlab('Dispersion of private values') + ylab('%') + geom_vline(xintercept = 1,linetype = 'dashed')
  ggsave('../out/figures/A12C.png',height=hh,width=ww,units = 'in',dpi = 300)
  
  
  
  
  # 1.a Extreme liquidity outcome
  liquidity.cs    <- runComparativeStatic(par = 'lambda',lower = .05, upper = 100, length = 3,calibrated.model)
  liquidity.cs    <- liquidity.cs[scale > .5]
  liquidity.cs$Scenario <- factor(x=c('Normal','Hot'),levels = c('Normal','Hot'))
  ggplot(liquidity.cs) + geom_bar(aes(x=Scenario,y=iBuyer_share),fill='black',color='black',width=.5,stat='identity') + theme_classic() + ylab(NULL) + xlab(NULL) + scale_y_continuous(labels = function(x) {paste0(100*x,'%')}) + geom_hline(yintercept =0)
  ggsave('../out/figures/A10A.png',height=hh,width=ww,units='in',dpi = 300)
  
  ggplot(liquidity.cs) + geom_bar(aes(x=Scenario,y=share_b_given_p),fill='black',color='black',width=.5,stat='identity') + theme_classic() + ylab(NULL) + xlab(NULL) + scale_y_continuous(labels = function(x) {paste0(100*x,'%')}) + geom_hline(yintercept =0)
  ggsave('../out/figures/A10B.png',height=hh,width=ww,units='in',dpi = 300)
  
  
  # Gross spread and fixed operating costs (Appendix Figure A.13)
  fixed_operating_cost = 0.01 * max(preference.cs$house_price_ib) 
  ggplot(preference.cs[iBuyer_share > 0.0175]) + 
    geom_line(aes(x=iBuyer_share,y=(profitability * house_price_ib * iBuyer_share * (365/tom_ib) - fixed_operating_cost)/(house_price_ib * iBuyer_share * (365/tom_ib))   )) +
    theme_classic() + xlab(NULL) + ylab(NULL) + geom_hline(yintercept = 0,color='gray') + scale_y_continuous(labels = function(x) {paste0(100*x,'%')},breaks = seq(-0.08,0.05,by=0.01))+ scale_x_continuous(labels = function(x) {paste0(100*x,'%')}) +
    coord_cartesian(ylim = c(-0.08,.05)) 
  ggsave('../out/figures/A13.png',height=4,width=6,units='in',dpi = 300)
  
}



getEmpiricalMoments_Pre_iBuyer <- function(X=1) {
  
  moments = rbind(data.table(name = 'tom.hh'             , value = 91,    weight = X * .1),
                  data.table(name = 'imp.price.delta.pct', value = 0.0565, weight =X *  260/100), 
                  data.table(name = 'imp.list.delta.days', value = 56.499,    weight = X / 100),
                  data.table(name = 'price_listed'       , value = 265,   weight = X * 1),
                  data.table(name = 'price_sd_xs'      , value =   0.180718, weight = X * 1), 
                  data.table(name = 'price_sd_ts'      , value =   0.09887161, weight = X * 1)) 
  
  
}


getEmpiricalMoments_Post_iBuyer <- function(X=1) {
  
  old.moments <- getEmpiricalMoments_Pre_iBuyer(X)
  
  new.moments = rbind(data.table(name = 'p.reduce.prices'    , value =   0.56197075,  weight = X * 1/0.56197075),      # From Listings Data
                      data.table(name = 'tom.ib.delta'       , value =  6          ,  weight = X * 1/6),               # From Listings Data
                      data.table(name = 'iBuyer.share'       , value =   0.04880825,  weight = X * 15/0.04880825),     # Market share as of 2018
                      data.table(name = 'buy.discount'       , value =  -.0264038,       weight = X * 5/.0264038),     # Table 2B Regression as of 2018
                      data.table(name = 'sell.premium'       , value =   0.0080326, weight = X * 0/0.016036897),       # Table 2B Regression as of 2018
                      data.table(name = 'diB.dLiq'           , value =   0.033434    , weight = 2/0.033434),           # From "pricing erros liquidity 4 struct march 2022" q6, 2018 only
                      data.table(name = 'diB.dErr'           , value =  -0.045495   , weight = X * 1/0.045495),        # From "pricing erros liquidity 4 struct march 2022" q3,  2018 only
                      data.table(name = 'DIST'               , value = 0            , weight = 999999999),             # Make sure it's solved
                      data.table(name = 'share.unprofitable' , value = 0.1278254          , weight =  10/0.1278254),       # From pnl_analysis_nov_2024 (fraction with pnl < 0), 
                      data.table(name = 'ibuyer.spread'      , value = 0.0344364   , weight = 500/0.056))              # Difference between buy discount and sale premium (target this, not sell premium)
  
  
  
  moments <- rbind(old.moments,new.moments)
  
  
  
}




setPars <- function(modifications = NULL) {
  
  
  
  
  # Deterministic parameters
  rho          = 0.050     # discount rate (External 1) (ASSUMED) # Bayer WP 2020 
  u_matched    =  15       # Normalization
  u_pat_delta  = 2         # How much worse unmatched is
  u_imp_delta  = 2         # How much EVEN worse unmatched impatient is
  u_ib_delta   = 6         # How much worse iBuyer is
  lambda       = 150       # Market liquidity (Cal 4) ~ Time on market
  lambda_ib    =  1        # iBuyers' relative skill in matching
  mu      = 9.1/60         # Separation rate (External 2) . https://www.census.gov/topics/population/migration/guidance/calculating-migration-expectancy.html . Move 9.1 times in 60 years --> rate = 9.1/60 
  M       = 1              # Population size (Fixed)
  k       = 0# .5772          # Negotiation cost (Fixed)
  ch      = 5              # Cost if roof caves in (Cal 5) ~ Change in iBuyer listing price 
  nu      = 1              # Probability that the roof caves in ~ (Cal 6) p(iBuyer updtes price)
  sig_trans = .14
  p.imp   = 0.5
  ibu        = 1           # HH hedonic utility of using ibuyer (TEST)
  ibuyer.overhead = 0.00
  
  # Volatility for preference shocks
  sig_m  = 1.0            # Volatility on liking the house. (Cal 7) ???
  sig_i  = 1.0            # Volatility on liking to sell to an iBuyer (Cal 8) ???
  
  # Adverse selection parameters
  # Shock parameters
  p_switch_types  = 0.25      # Probability that house is a bad type
  signal_noise    = 0.01      # Noise term (Cal 9) Relationship between iBuyer and pricing
  u_delta_bad     = 2         # Change in flow utility for (human) homeowner if its a bad house 
  dispersion      = 0
  
  # Delay in iBuyer closing
  delay = 15/365
  
  if(!is.null(modifications)) {
    for(ff in 1:nrow(modifications)) {
      eval(parse(text=paste0(modifications$variable[ff],'=',modifications$value[ff])))
      
    }
  }
  
  eps = 1e-30
  
  
  u_matched_good = u_matched + dispersion
  u_matched_bad  = u_matched     - u_delta_bad - dispersion
  u_unm_pat      = u_matched     - u_pat_delta # WHY DELTA BAD HERE
  u_unm_imp      = u_matched     - u_pat_delta - u_imp_delta  # WHY DELTA BAD HERE
  u_ibuyer       = u_matched - u_ib_delta
  u_ibuyer_shock = u_matched - u_ib_delta  - ch
  
  fb =  function(b,s) { lambda * s }           # Matching function for individual buyer
  fs =  function(b,s) { lambda * b }           # Matching function for individual seller
  f  =  function(b,s) { lambda * b * s }       # Aggregate matching function
  
  
  
  pars <- list(ibu = ibu,u_ibuyer_shock = u_ibuyer_shock,signal_noise = signal_noise,ibuyer.overhead = ibuyer.overhead,sig_m = sig_m, sig_i = sig_i, rho = rho,p.imp = p.imp,sig_trans = sig_trans, delay = delay,u_unm_pat = u_unm_pat,u_unm_imp =u_unm_imp, u_matched_good = u_matched_good,u_matched_bad = u_matched_bad,p_switch_types = p_switch_types, u_ibuyer = u_ibuyer, lambda = lambda, lambda_ib = lambda_ib, mu = mu, M = M,k = k, fb = fb, fs = fs, f = f,ch = ch,nu = nu,eps = eps)
  
  return(pars)
}







runComparativeStatic <- function(par = 'lambda',lower = .1,upper = 10,length = 100,calibrated.model) {
  
  
  
  par.space = c(seq(lower,1,length.out=length/2),seq(1,upper,length.out = length/2))
  baseline.parameters = calibrated.model$par.dt
  baseline.parameters <- rbind(baseline.parameters,data.table(variable = 'delay',value = 15/365))
  baseline.parameters <- rbind(baseline.parameters,data.table(variable = 'dispersion',value = 0))
  

  
  baseline.no.ibuyer.guess = calibrated.model$no.ibuyer.guess
  baseline.ibuyer.guess = calibrated.model$ibuyer.guess
  no.ibuyer.guess <- calibrated.model$no.ibuyer.guess
  ibuyer.guess <- calibrated.model$ibuyer.guess
  
  
  if('dispersion' %in% par) {
    up.range = c(0,par.space[par.space > 0]) 
  } else {
    up.range = c(1,par.space[par.space > 1])
  }
  temp.guess <<- ibuyer.guess
  for(pp in 1:length(up.range)) {
    
    # First go up from solved, then go down
    current.parameters <- baseline.parameters
    current.parameters[variable %in% par]$value <-   current.parameters[variable %in% par]$value * up.range[pp]
    if('dispersion' %in% par) {
      current.parameters[variable == 'dispersion']$value <-   up.range[pp] 
    }
    print(current.parameters[variable %in% par]$value)
    pars <- setPars(current.parameters)
    solved <- solveModelRevision(pars,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = temp.guess,pre_ibuyer_only = F) 
    dist <- solved$ibuyer$stats$DIST
    if(dist < 1e-9) {
      temp.guess <<- solved$ibuyer$solution
    }
    
    share_b_given_g = sum(solved$ibuyer$stats$ib_prod_share[past == 'g' & signal == 'p' & type == 'b']$share)
    share_b         = sum(solved$ibuyer$stats$ib_prod_share[type == 'b']$share)
    share_b_given_p = sum(solved$ibuyer$stats$ib_prod_share[signal == 'p' & type == 'b']$share) / sum(solved$ibuyer$stats$ib_prod_share[signal == 'p']$share)
    profitability   = sum(solved$ibuyer$stats$ib_prod_share$share * solved$ibuyer$stats$ib_prod_share$profitability)/solved$no.ibuyer$stats$price_hh_mean
    share_pat       = sum(solved$ibuyer$stats$ib_prod_share[patient  == 'p']$share)
    share_imp       = sum(solved$ibuyer$stats$ib_prod_share[patient  == 'i']$share)
    tom             = solved$ibuyer$stats$tom_ib_mean
    hp_ib           = solved$ibuyer$stats$price_hh_mean
    hp_nib          = solved$no.ibuyer$stats$price_hh_mean
    welfare_ib      = solved$no.ibuyer$stats$vhh
    welfare_nib     = solved$ibuyer$stats$vhh
    
    
   
    
    
    
    temp.result <- data.table(scale = up.range[pp],
                              param.value = current.parameters[variable %in% par]$value[1],
                              iBuyer_share = solved$ibuyer$stats$ib_share,
                              mispriced_share = share_b_given_g,
                              share_b_given_p = share_b_given_p,
                              price.vol = solved$no.ibuyer$stats$log_price_hh_std_xs,
                              tom = solved$no.ibuyer$stats$tom_hh_mean,
                              profitability = profitability,
                              share_b = share_b,
                              share_pat = share_pat,
                              share_imp = share_imp,
                              tom_ib = tom,
                              house_price_ib = hp_ib,
                              house_price_nib = hp_nib,
                              welfare_ib = welfare_ib,
                              welfare_nib = welfare_nib,
                              v_ib = solved$ibuyer$stats$v_ib_mean,
                              dist = dist)
    
    if(pp == 1) {
      result = temp.result
    } else {
      result = rbind(result,temp.result)
    }
  }
  
  if('dispersion' %in% par) {
    down.range =  c(0,par.space[par.space < 0]) 
    down.range <- down.range[order(-down.range)]
    
  } else {
    down.range <- par.space[par.space < 1]
    down.range <- down.range[order(-down.range)]
    
  }
  temp.guess <<- calibrated.model$ibuyer.guess
  for(pp in 1:length(down.range)) {
    
    # First go up from solved, then go down
    current.parameters <- baseline.parameters
    current.parameters[variable %in% par]$value <-   current.parameters[variable %in% par]$value * down.range[pp]
    if('dispersion' %in% par) {
      current.parameters[variable == 'dispersion']$value <-   down.range[pp] 
    }
    print(current.parameters[variable %in% par]$value)
    pars <- setPars(current.parameters)
    solved <- solveModelRevision(pars,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = temp.guess,pre_ibuyer_only = F) 
    dist <- solved$ibuyer$stats$DIST
    if(dist < 1e-15) {
      temp.guess <<- solved$ibuyer$solution
    }
    
    share_b_given_g = sum(solved$ibuyer$stats$ib_prod_share[past == 'g' & signal == 'p' & type == 'b']$share)
    share_b         = sum(solved$ibuyer$stats$ib_prod_share[type == 'b']$share)
    profitability   = sum(solved$ibuyer$stats$ib_prod_share$share * solved$ibuyer$stats$ib_prod_share$profitability)/solved$no.ibuyer$stats$price_hh_mean
    share_b_given_p = sum(solved$ibuyer$stats$ib_prod_share[signal == 'p' & type == 'b']$share) / sum(solved$ibuyer$stats$ib_prod_share[signal == 'p']$share)
    share_pat       = sum(solved$ibuyer$stats$ib_prod_share[patient  == 'p']$share)
    share_imp       = sum(solved$ibuyer$stats$ib_prod_share[patient  == 'i']$share)
    tom             = solved$ibuyer$stats$tom_ib_mean
    hp_ib           = solved$ibuyer$stats$price_hh_mean
    hp_nib          = solved$no.ibuyer$stats$price_hh_mean
    welfare_ib      = solved$no.ibuyer$stats$vhh
    welfare_nib     = solved$ibuyer$stats$vhh
    
    temp.result <- data.table(scale = down.range[pp],
                              param.value = current.parameters[variable %in% par]$value[1],
                              iBuyer_share = solved$ibuyer$stats$ib_share,
                              mispriced_share = share_b_given_g,
                              share_b_given_p = share_b_given_p,
                              tom = solved$no.ibuyer$stats$tom_hh_mean,
                              price.vol = solved$no.ibuyer$stats$log_price_hh_std_xs,
                              profitability = profitability,
                              share_b = share_b,
                              share_pat = share_pat,
                              share_imp = share_imp,
                              tom_ib = tom,
                              house_price_ib = hp_ib,
                              house_price_nib = hp_nib,
                              welfare_ib = welfare_ib,
                              welfare_nib = welfare_nib,
                              v_ib = solved$ibuyer$stats$v_ib_mean,
                              dist = dist)
    
    result = rbind(result,temp.result)
  }
  
  return(result)
  
}



runComparativeStaticDT <- function(dt,calibrated.model) {
  
  baseline.parameters = calibrated.model$par.dt
  baseline.parameters <- rbind(baseline.parameters,data.table(variable = 'delay',value = 15/365))
  baseline.parameters <- baseline.parameters[!(variable %in% dt$variable)]
  new.parameters <- rbind(baseline.parameters,dt)
  
  baseline.no.ibuyer.guess = calibrated.model$no.ibuyer.guess
  baseline.ibuyer.guess = calibrated.model$ibuyer.guess
  no.ibuyer.guess <- calibrated.model$no.ibuyer.guess
  ibuyer.guess <- calibrated.model$ibuyer.guess
  
  pars <- setPars(new.parameters)
  solved <- solveModelRevision(pars,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = ibuyer.guess,pre_ibuyer_only = F) 
  
  
    share_b_given_g = sum(solved$ibuyer$stats$ib_prod_share[past == 'g' & signal == 'p' & type == 'b']$share)
    share_b         = sum(solved$ibuyer$stats$ib_prod_share[type == 'b']$share)
    share_b_given_p = sum(solved$ibuyer$stats$ib_prod_share[signal == 'p' & type == 'b']$share) / sum(solved$ibuyer$stats$ib_prod_share[signal == 'p']$share)
    profitability   = sum(solved$ibuyer$stats$ib_prod_share$share * solved$ibuyer$stats$ib_prod_share$profitability)/solved$no.ibuyer$stats$price_hh_mean
    share_pat       = sum(solved$ibuyer$stats$ib_prod_share[patient  == 'p']$share)
    share_imp       = sum(solved$ibuyer$stats$ib_prod_share[patient  == 'i']$share)
    
    temp.result <- data.table(
                              iBuyer_share = solved$ibuyer$stats$ib_share,
                              mispriced_share = share_b_given_g,
                              share_b_given_p = share_b_given_p,
                              price.vol = solved$no.ibuyer$stats$log_price_hh_std_xs,
                              tom = solved$no.ibuyer$stats$tom_hh_mean,
                              profitability = profitability,
                              share_b = share_b,
                              share_pat = share_pat,
                              share_imp = share_imp,
                              dist = dist)
    
   
  
  return(temp.result)
  
}

getIBShareDerivatives <- function(pars.dt,no.ibuyer.guess = NULL,ibuyer.guess = NULL) {
  
  ### 1. Change noise
  base.pars <- setPars(pars.dt)
  
  
  # DELTA 
  delta = .01
  
  m0 <- solveModelRevision(base.pars,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = ibuyer.guess,pre_ibuyer_only = F)
  
  pars.dt.1 <- pars.dt
  pars.dt.1[variable == 'signal_noise']$value <- pars.dt[variable == 'signal_noise']$value + delta 
  pars1 <- setPars(pars.dt.1)
  
  m1 <- solveModelRevision(pars1,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = m0$ibuyer$solution,pre_ibuyer_only = F)
  
  s0 = m0$ibuyer$stats$ib_share
  s1 = m1$ibuyer$stats$ib_share
  
  sd0 = m0$no.ibuyer$stats$log_price_hh_std_xs + pars.dt[variable == 'signal_noise']$value
  sd1 = m1$no.ibuyer$stats$log_price_hh_std_xs + pars.dt.1[variable == 'signal_noise']$value
  
  # Derivative
  ds_dn = (s1-s0)/(sd1-sd0)
  
  
  pars.dt.2 <- pars.dt
  pars.dt.2[variable == 'lambda']$value <- pars.dt[variable == 'lambda']$value + delta
  pars.dt.2[variable == 'p.imp']$value <- pars.dt[variable == 'p.imp']$value + .1
  
  pars2 <- setPars(pars.dt.2)
  
  
  m2 <- solveModelRevision(pars2,no.ibuyer.guess = no.ibuyer.guess,ibuyer.guess = m0$ibuyer$solution,pre_ibuyer_only = F)
  
  s0 = m0$ibuyer$stats$ib_share
  s2 = m2$ibuyer$stats$ib_share
  
  sd0 = m0$no.ibuyer$stats$p.sells.90.days
  sd2 = m2$no.ibuyer$stats$p.sells.90.days
  
  # Derivative
  ds_dv = (s2-s0)/(sd2-sd0)
  
  
  return(list(ds_dn = ds_dn,ds_dv = ds_dv))
  
}



solveModelRevision <- function(pars,no.ibuyer.guess = NULL,ibuyer.guess = NULL,pre_ibuyer_only = F) {
  
  
  # Solve the model when there's no iBuyer
  solveNoiBuyer <- function(x,return_stats = F) {
    p_hh_good_pat = x[1]
    p_hh_good_imp = x[2]
    p_hh_bad_pat  = x[3]
    p_hh_bad_imp  = x[4]
    vs_good_pat   = x[5] 
    vs_good_imp   = x[6]
    vs_bad_pat    = x[7]
    vs_bad_imp    = x[8]
    vb            = x[9]
    m_h_g         = x[10]
    m_h_b         = x[11]
    
    # Setup 
    sig_m         = pars$sig_m
    
    # Expected selling value functions before we know if you're impatient or not
    vs_good_expected = pars$p.imp * vs_good_imp + (1 - pars$p.imp) * vs_good_pat
    vs_bad_expected  = pars$p.imp * vs_bad_imp  + (1 - pars$p.imp) * vs_bad_pat
    
    # Directly calculate value of being well-matched for good and bad type houses
    vh_good = (pars$u_matched_good + pars$mu * ( (pars$p_switch_types) * vs_bad_expected  + (1 - pars$p_switch_types) * vs_good_expected ) ) / ( pars$rho + pars$mu)
    vh_bad  = (pars$u_matched_bad  + pars$mu * ( (pars$p_switch_types) * vs_good_expected + (1 - pars$p_switch_types) * vs_bad_expected  ) ) / ( pars$rho + pars$mu)
    
    # Directly calculate values when you buy from various types
    v_buy_good_from_pat = vh_good - vb - p_hh_good_pat
    v_buy_good_from_imp = vh_good - vb - p_hh_good_imp
    v_buy_bad_from_pat  = vh_bad  - vb - p_hh_bad_pat
    v_buy_bad_from_imp  = vh_bad  - vb - p_hh_bad_imp
    
    # Calculate probabilities that you buy from a given offer
    pi_buy_good_from_pat   = 1   / ( exp(-v_buy_good_from_pat / sig_m) + 1 ) 
    pi_buy_good_from_imp   = 1   / ( exp(-v_buy_good_from_imp / sig_m) + 1 ) 
    pi_buy_bad_from_pat    = 1   / ( exp(-v_buy_bad_from_pat  / sig_m) + 1 ) 
    pi_buy_bad_from_imp    = 1   / ( exp(-v_buy_bad_from_imp  / sig_m) + 1 ) 
    
    # Calculate laws of motion
    m_b    = (pars$M - m_h_g - m_h_b)/2
    m_s_pg = pars$mu * (1 - pars$p.imp) * (m_h_g * (1 - pars$p_switch_types) + m_h_b * (    pars$p_switch_types)  ) / (pars$lambda * m_b * pi_buy_good_from_pat)
    m_s_pb = pars$mu * (1 - pars$p.imp) * (m_h_g * (    pars$p_switch_types) + m_h_b * (1 - pars$p_switch_types)  ) / (pars$lambda * m_b * pi_buy_bad_from_pat )
    m_s_ig = pars$mu * (    pars$p.imp) * (m_h_g * (1 - pars$p_switch_types) + m_h_b * (    pars$p_switch_types)  ) / (pars$lambda * m_b * pi_buy_good_from_imp)
    m_s_ib = pars$mu * (    pars$p.imp) * (m_h_g * (    pars$p_switch_types) + m_h_b * (1 - pars$p_switch_types)  ) / (pars$lambda * m_b * pi_buy_bad_from_imp )
    
    # Expected value when you meet with a seller (even if you don't end up buying from her) # Why is k multiplied by sig_m? Probably should be outside.
    Eu_buy_from_good_imp    = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_good_from_imp))  # Become happy, stop being buyer, pay impatient hh price
    Eu_buy_from_good_pat    = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_good_from_pat))  # Become happy, stop being buyer, pay patient hh price
    Eu_buy_from_bad_imp     = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_bad_from_imp))  # Become happy, stop being buyer, pay impatient hh price
    Eu_buy_from_bad_pat     = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_bad_from_pat))  # Become happy, stop being buyer, pay patient hh price
    
    # Conditions to zero out
    
    # Value functions
    eqn_vs_good_pat       = -vs_good_pat             + (pars$u_unm_pat  +  pars$fs(b = m_b,s = 1)    * (pi_buy_good_from_pat * (p_hh_good_pat + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_good_from_pat) # You may find a guy to sell to :)
    eqn_vs_good_imp       = -vs_good_imp             + (pars$u_unm_imp  +  pars$fs(b = m_b,s = 1)    * (pi_buy_good_from_imp * (p_hh_good_imp + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_good_from_imp) # You may find a guy to sell to :)
    eqn_vs_bad_pat        = -vs_bad_pat              + (pars$u_unm_pat  +  pars$fs(b = m_b,s = 1)    * (pi_buy_bad_from_pat  * (p_hh_bad_pat  + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_bad_from_pat) # You may find a guy to sell to :)
    eqn_vs_bad_imp        = -vs_bad_imp              + (pars$u_unm_imp  +  pars$fs(b = m_b,s = 1)    * (pi_buy_bad_from_imp  * (p_hh_bad_imp  + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_bad_from_imp) # You may find a guy to sell to :)
    eqn_vb                = -vb             
                            + (pars$u_unm_pat  +  
                                 pars$fb(b = 1,s = m_s_pg) * Eu_buy_from_good_pat  + 
                                 pars$fb(b = 1,s = m_s_pb) * Eu_buy_from_bad_pat   + 
                                 pars$fb(b = 1,s = m_s_ig) * Eu_buy_from_good_imp  + 
                                 pars$fb(b = 1,s = m_s_ib) * Eu_buy_from_bad_imp  ) / pars$rho  
    
    # First-order conditions for pricing
    eqn_foc_s_good_pat    = - 1 / sig_m * (1 - pi_buy_good_from_pat)   * (vb + p_hh_good_pat - vs_good_pat)   + 1 # Homeowner's FOC when listing
    eqn_foc_s_good_imp    = - 1 / sig_m * (1 - pi_buy_good_from_imp)   * (vb + p_hh_good_imp - vs_good_imp)   + 1 # Homeowner's FOC when listing
    eqn_foc_s_bad_pat     = - 1 / sig_m * (1 - pi_buy_bad_from_pat)    * (vb + p_hh_bad_pat  - vs_bad_pat )   + 1 # Homeowner's FOC when listing
    eqn_foc_s_bad_imp     = - 1 / sig_m * (1 - pi_buy_bad_from_imp)    * (vb + p_hh_bad_imp  - vs_bad_imp )   + 1 # Homeowner's FOC when listing
    
    # Buyers = sellers
    eqn_buyers_equal_sellers = m_b - m_s_pg - m_s_pb - m_s_ig - m_s_ib
    
    # Equal number of good and bad houses
    eqn_good_bad_houses = (m_h_g + m_s_pg + m_s_ig) - (m_h_b + m_s_pb + m_s_ib)
    if(!return_stats) { 
      return(c(eqn_vs_good_pat,
               eqn_vs_good_imp,
               eqn_vs_bad_pat,
               eqn_vs_bad_imp,
               eqn_vb,
               eqn_foc_s_good_pat,
               eqn_foc_s_good_imp,
               eqn_foc_s_bad_pat,
               eqn_foc_s_bad_imp,
               eqn_buyers_equal_sellers,
               eqn_good_bad_houses))
    } else {
      
      stats = list()
      
      # Value functions
      stats$vs_good_pat = vs_good_pat
      stats$vs_good_imp = vs_good_imp
      stats$vs_bad_pat  = vs_bad_pat
      stats$vs_bad_imp  = vs_bad_imp
      stats$vb          = vb
      stats$vh_good     = vh_good
      stats$vh_bad      = vh_bad
      stats$vhh         = m_h_g * vh_good + m_h_b * vh_bad + m_b * vb + m_s_pg * vs_good_pat + m_s_pb * vs_bad_pat + m_s_ig * vs_good_imp + m_s_ib * vs_bad_imp
      
      # Prices
      flow_rate_hh_2_hh_g_p = pars$lambda * m_b * m_s_pg * pi_buy_good_from_pat
      flow_rate_hh_2_hh_g_i = pars$lambda * m_b * m_s_ig * pi_buy_good_from_imp
      flow_rate_hh_2_hh_b_p = pars$lambda * m_b * m_s_pb * pi_buy_bad_from_pat 
      flow_rate_hh_2_hh_b_i = pars$lambda * m_b * m_s_ib * pi_buy_bad_from_imp 
      hh_2_matched_flow_rate = flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_g_i + flow_rate_hh_2_hh_b_p + flow_rate_hh_2_hh_b_i
      
      
      stats$price_hh_good_patient   = p_hh_good_pat
      stats$price_hh_good_impatient = p_hh_good_imp
      stats$price_hh_bad_patient    = p_hh_bad_pat
      stats$price_hh_bad_impatient  = p_hh_bad_imp
      stats$price_hh_mean           = 1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_p * p_hh_good_pat   + flow_rate_hh_2_hh_g_i * p_hh_good_imp   + flow_rate_hh_2_hh_b_p * p_hh_bad_pat   + flow_rate_hh_2_hh_b_i * p_hh_bad_imp  )
      log_price_hh_mean             = 1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_p * log(p_hh_good_pat)   + flow_rate_hh_2_hh_g_i * log(p_hh_good_imp)   + flow_rate_hh_2_hh_b_p * log(p_hh_bad_pat)   + flow_rate_hh_2_hh_b_i * log(p_hh_bad_imp)  )
      v_log_price_hh                = 1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_p * log(p_hh_good_pat)^2 + flow_rate_hh_2_hh_g_i * log(p_hh_good_imp)^2 + flow_rate_hh_2_hh_b_p * log(p_hh_bad_pat)^2 + flow_rate_hh_2_hh_b_i * log(p_hh_bad_imp)^2) 
      stats$log_price_hh_std_xs     = sqrt( v_log_price_hh - log_price_hh_mean^2  )#     + (3.14/sqrt(6)*pars$sig_m/stats$price_hh_mean)^2)
      stats$price_hh_mean_patient   = 1 / (flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_b_p) * (flow_rate_hh_2_hh_g_p * p_hh_good_pat   + flow_rate_hh_2_hh_b_p * p_hh_bad_pat  )
      stats$price_hh_mean_impatient = 1 / (flow_rate_hh_2_hh_g_i + flow_rate_hh_2_hh_b_i) * (flow_rate_hh_2_hh_g_i * p_hh_good_imp   + flow_rate_hh_2_hh_b_i * p_hh_bad_imp  )
      
      
      
      # Calculate time-series standard deviation in logs
      
      # 1. Get var | last time
      Ep_gp = (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_pat)   + (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_pat)   + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_imp)   + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_imp)
      Ep_gi = (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_pat)   + (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_pat)   + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_imp)   + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_imp)
      Ep_bp = (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_pat)   + (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_pat)   + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_imp)   + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_imp)
      Ep_bi = (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_pat)   + (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_pat)   + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_imp)   + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_imp)
      
      Vp_gp = (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_pat)^2 + (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_pat)^2 + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_imp)^2 + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_imp)^2
      Vp_gi = (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_pat)^2 + (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_pat)^2 + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_good_imp)^2 + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_bad_imp)^2
      Vp_bp = (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_pat)^2 + (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_pat)^2 + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_imp)^2 + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_imp)^2
      Vp_bi = (1 - pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_pat)^2 + (1 - pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_pat)^2 + (    pars$p.imp) * (    pars$p_switch_types) * log(p_hh_good_imp)^2 + (    pars$p.imp) * (1 - pars$p_switch_types) * log(p_hh_bad_imp)^2
      
      stats$log_price_hh_std_ts  = sqrt(1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_p * (Vp_gp - Ep_gp^2) + flow_rate_hh_2_hh_g_i * (Vp_gi - Ep_gp^2)
                                                                    + flow_rate_hh_2_hh_b_p * (Vp_bp - Ep_bp^2) + flow_rate_hh_2_hh_b_i * (Vp_bi - Ep_bi^2)) )# + (3.14/sqrt(6)*pars$sig_m/stats$price_hh_mean)^2)
      
      
      # Time on market
      stats$tom_hh_good_patient     = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_good_from_pat) * 365
      stats$tom_hh_good_impatient   = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_good_from_imp) * 365 
      stats$tom_hh_bad_patient      = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_bad_from_pat ) * 365 
      stats$tom_hh_bad_impatient    = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_bad_from_imp ) * 365 
      stats$tom_hh_mean             = 1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_i * stats$tom_hh_good_impatient + flow_rate_hh_2_hh_b_i * stats$tom_hh_bad_impatient + flow_rate_hh_2_hh_g_p * stats$tom_hh_good_patient + flow_rate_hh_2_hh_b_p * stats$tom_hh_bad_patient) 
      stats$tom_hh_patient_mean     = 1 / (flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_b_p) * (flow_rate_hh_2_hh_g_p * stats$tom_hh_good_patient   + flow_rate_hh_2_hh_b_p * stats$tom_hh_bad_patient  )
      stats$tom_hh_impatient_mean   = 1 / (flow_rate_hh_2_hh_g_i + flow_rate_hh_2_hh_b_i) * (flow_rate_hh_2_hh_g_i * stats$tom_hh_good_impatient + flow_rate_hh_2_hh_b_i * stats$tom_hh_bad_impatient)
      
      
      # Search time (how long buyer spends looking)
      stats$search_time             = 1 / (pars$fb(b = 1,s = m_s_ig)  * pi_buy_good_from_imp   + 
                                           pars$fb(b = 1,s = m_s_pg)  * pi_buy_good_from_pat   + 
                                           pars$fb(b = 1,s = m_s_ib)  * pi_buy_bad_from_imp    + 
                                           pars$fb(b = 1,s = m_s_pb)  * pi_buy_bad_from_pat)   * 365 
      
      
      stats$p.sells.90.days = 1 - exp(-0.25 * 1 / (stats$tom_hh_mean / 365) )
      
      
      
      return(stats)
      
    }
  }
  
  
  
  # 1. Help out by bootstrapping a no-ibuyer solution and using that as initial values for the iBuyer problem.
  if(!is.null(no.ibuyer.guess)) {
    guess.no.ibuyer <- no.ibuyer.guess
  } else {
    guess.no.ibuyer = c(pars$u_matched_good / pars$rho , pars$u_matched_good / pars$rho - 1, pars$u_matched_bad / pars$rho, pars$u_matched_bad / pars$rho - 1,
                        pars$u_matched_good / pars$rho , pars$u_matched_good / pars$rho - 1, pars$u_matched_bad / pars$rho, pars$u_matched_bad / pars$rho - 1,
                        0,.45,.45)
  }
  
  no.ibuyer.soln <- nleqslv(x = guess.no.ibuyer,fn = solveNoiBuyer,method = 'Broyden',xscalm = 'fixed',control = list(maxit = 50000,xtol = 1e-10,ftol = 1e-10,stepmax = .1,allowSingular = T))
  no.ibuyer.solution = no.ibuyer.soln$x
  no.ibuyer.stats <- solveNoiBuyer(no.ibuyer.solution,return_stats = T)
  
  
  if(!pre_ibuyer_only) {
    
    ## Solve it with iBuyers
    
    # Actual equilibrium to solve
    toSolve <- function(x,return_stats = F) {
      
      
      # Household listing prices (x4: good/bad x patient/impatient)
      p_hh_good_pat = x[1]
      p_hh_good_imp = x[2]
      p_hh_bad_pat  = x[3]
      p_hh_bad_imp  = x[4]
      
      # iBuyer prices when buying (x4: positive/negative signal x good/bad last time)
      p_ib_b_p_g    = x[5]# p for "positive" signal, n for negative signal; g for good last time, b for bad last time
      p_ib_b_p_b    = x[6]
      p_ib_b_n_g    = x[7]
      p_ib_b_n_b    = x[8]
      
      # iBuyer prices when selling (x4: good/bad house x no-shock/shock)
      p_ib_s_g_e    = x[9] # g/b for good/bad (this time); e for "early" (pre-roof leak), l for "late" (post-roof leak)
      p_ib_s_g_l    = x[10]
      p_ib_s_b_e    = x[11]
      p_ib_s_b_l    = x[12]
      
      # Value functions (x5: 4x seller + 1x buyer)
      vs_good_pat   = x[13] 
      vs_good_imp   = x[14]
      vs_bad_pat    = x[15]
      vs_bad_imp    = x[16]
      
      # Masses
      m_h_g         = x[17]
      m_h_b         = x[18]
      m_ib_g_e      = x[19] 
      m_ib_g_l      = x[20] 
      m_ib_b_e      = x[21]
      m_ib_b_l      = x[22]
      
      # Implied by eqm.
      vb = 0
      
      # Some parameters that show up a lot
      sig_m = pars$sig_m
      sig_i = pars$sig_i
      p_switch_types = pars$p_switch_types
      signal_noise   = pars$signal_noise
      p.imp          = pars$p.imp
      ddf = exp(-pars$delay * pars$rho)
      
      # Directly calculate some value functions
      
      # Values of selling to iBuyers---depends on:
      # -- was it good or bad last time (x2) (iBuyer knows)
      # -- did you get a pos/neg signal (x2) (iBuyer knows)
      # -- are you patient or impatient (x2) (iBuyer doesn't know)
      # -- is it good or bad now (x2)        (iBuyer doesn't know) (but this doesn't matter for value of selling to iBuyer)
      #
      # -----> 8 values
      
      # _g_p_p: (house used to be good; you got a positive signal; you're patient)
      v_sellToiBuyer_g_p_p = (vb + p_ib_b_p_g) * ddf + pars$u_unm_pat / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu
      v_sellToiBuyer_g_p_i = (vb + p_ib_b_p_g) * ddf + pars$u_unm_imp / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu
      
      v_sellToiBuyer_g_n_p = (vb + p_ib_b_n_g) * ddf + pars$u_unm_pat / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu 
      v_sellToiBuyer_g_n_i = (vb + p_ib_b_n_g) * ddf + pars$u_unm_imp / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu 
        
      v_sellToiBuyer_b_p_p = (vb + p_ib_b_p_b) * ddf + pars$u_unm_pat / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu
      v_sellToiBuyer_b_p_i = (vb + p_ib_b_p_b) * ddf + pars$u_unm_imp / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu
        
      v_sellToiBuyer_b_n_p = (vb + p_ib_b_n_b) * ddf + pars$u_unm_pat / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu
      v_sellToiBuyer_b_n_i = (vb + p_ib_b_n_b) * ddf + pars$u_unm_imp / pars$rho * (1 - exp(-pars$rho * pars$delay)) + pars$ibu
      
      
      # What is the matched-type's value when she gets a shock?
      # If your old house was a good house
      Eu_receive_shock_g  = sig_i * .5772 +           # Option value from T1EV shock
                            (1-p.imp) * (             # You're patient
                              (1 - p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_g_p_p) + exp(vs_good_pat)) +  # Stay good type and get positive signal
                              (1 - p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_g_n_p) + exp(vs_good_pat)) +  # Stay good type and get negative signal
                              (    p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_g_p_p) + exp(vs_bad_pat )) +  # Goto bad  type and get positive signal
                              (    p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_g_n_p) + exp(vs_bad_pat ))    # Goto bad  type and get negative signal
                            ) +
                            (  p.imp) * (             # You're impatient
                              (1 - p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_g_p_i) + exp(vs_good_imp)) +  # Stay good type and get positive signal
                              (1 - p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_g_n_i) + exp(vs_good_imp)) +  # Stay good type and get negative signal
                              (    p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_g_p_i) + exp(vs_bad_imp))  +  # Goto bad  type and get positive signal
                              (    p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_g_n_i) + exp(vs_bad_imp))    # Goto bad  type and get negative signal
                            )
      
      # If your old house was a bad house
      Eu_receive_shock_b  = sig_i * .5772 +           # Option value from T1EV shock
                            (1-p.imp) * (             # You're patient
                              (    p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_b_p_p) + exp(vs_good_pat)) +  # Goto good type and get positive signal
                              (    p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_b_n_p) + exp(vs_good_pat)) +  # Goto good type and get negative signal
                              (1 - p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_b_p_p) + exp(vs_bad_pat )) +  # Stay bad  type and get positive signal
                              (1 - p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_b_n_p) + exp(vs_bad_pat ))    # Stay bad  type and get negative signal
                            ) +
                            (  p.imp) * (             # You're impatient
                              (    p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_b_p_i) + exp(vs_good_imp)) +  # Goto good type and get positive signal
                              (    p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_b_n_i) + exp(vs_good_imp)) +  # Goto good type and get negative signal
                              (1 - p_switch_types) * (    signal_noise) * log(exp(v_sellToiBuyer_b_p_i) + exp(vs_bad_imp))  +  # Stay bad  type and get positive signal
                              (1 - p_switch_types) * (1 - signal_noise) * log(exp(v_sellToiBuyer_b_n_i) + exp(vs_bad_imp))     # Stay bad  type and get negative signal
                            )
      
      v_h_g             = (pars$u_matched_good  +  pars$mu * Eu_receive_shock_g) / (pars$rho + pars$mu) # You get shock with probability mu.
      v_h_b             = (pars$u_matched_bad   +  pars$mu * Eu_receive_shock_b) / (pars$rho + pars$mu) # You get shock with probability mu.
      
      
      # Value of transitioning from various sellers (impatient HH, patient HH, new iBuyers, old iBuyers)
      v_buy_from_good_pat = v_h_g - vb - p_hh_good_pat
      v_buy_from_good_imp = v_h_g - vb - p_hh_good_imp
      v_buy_from_bad_pat  = v_h_b - vb - p_hh_bad_pat
      v_buy_from_bad_imp  = v_h_b - vb - p_hh_bad_imp
        
      v_buy_from_good_ib_e = v_h_g - vb - p_ib_s_g_e
      v_buy_from_good_ib_l = v_h_g - vb - p_ib_s_g_l
      v_buy_from_bad_ib_e  = v_h_b - vb - p_ib_s_b_e
      v_buy_from_bad_ib_l  = v_h_b - vb - p_ib_s_b_l
      
      # Probabilities of buying from each type upon meeting
      pi_buy_from_good_pat = 1 / (exp(-v_buy_from_good_pat   / sig_m) + 1)
      pi_buy_from_good_imp = 1 / (exp(-v_buy_from_good_imp   / sig_m) + 1)
      pi_buy_from_bad_pat  = 1 / (exp(-v_buy_from_bad_pat    / sig_m) + 1)
      pi_buy_from_bad_imp  = 1 / (exp(-v_buy_from_bad_imp    / sig_m) + 1)
      
      pi_buy_from_good_ib_e = 1 / (exp(-v_buy_from_good_ib_e   / sig_m) + 1)
      pi_buy_from_good_ib_l = 1 / (exp(-v_buy_from_good_ib_l   / sig_m) + 1)
      pi_buy_from_bad_ib_e  = 1 / (exp(-v_buy_from_bad_ib_e    / sig_m) + 1)
      pi_buy_from_bad_ib_l  = 1 / (exp(-v_buy_from_bad_ib_l    / sig_m) + 1)
     
      # Probabilities of good/bad conditional on past type and signal
      pg_g_p = (1 - p_switch_types) * (1 - signal_noise) / ((1 - p_switch_types) * (1 - signal_noise) + (    p_switch_types) * (    signal_noise)  ) # Probability that it's good given you were good before and got the positive signal
      pg_g_n = (1 - p_switch_types) * (    signal_noise) / ((1 - p_switch_types) * (    signal_noise) + (    p_switch_types) * (1 - signal_noise)  ) # Probability that it's good given you were good before and got the negative signal
      pg_b_p = (    p_switch_types) * (1 - signal_noise) / ((    p_switch_types) * (1 - signal_noise) + (1 - p_switch_types) * (    signal_noise)  ) # Probability that it's good given you were good before and got the negative signal
      pg_b_n = (    p_switch_types) * (    signal_noise) / ((    p_switch_types) * (    signal_noise) + (1 - p_switch_types) * (1 - signal_noise)  ) # Probability that it's good given you were good before and got the negative signal
      
      # Households' decision to sell to iBuyer for (past type | signal | patient | current type).
      pi_sell_to_ibuyer_g_p_p_g = exp( (v_sellToiBuyer_g_p_p - vs_good_pat) / sig_i ) / ( exp( (v_sellToiBuyer_g_p_p - vs_good_pat) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_g_p_p_b = exp( (v_sellToiBuyer_g_p_p - vs_bad_pat ) / sig_i ) / ( exp( (v_sellToiBuyer_g_p_p - vs_bad_pat ) / sig_i ) + 1 ) 
      pi_sell_to_ibuyer_g_p_i_g = exp( (v_sellToiBuyer_g_p_i - vs_good_imp) / sig_i ) / ( exp( (v_sellToiBuyer_g_p_i - vs_good_imp) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_g_p_i_b = exp( (v_sellToiBuyer_g_p_i - vs_bad_imp ) / sig_i ) / ( exp( (v_sellToiBuyer_g_p_i - vs_bad_imp ) / sig_i ) + 1 )    
      
      pi_sell_to_ibuyer_g_n_p_g = exp( (v_sellToiBuyer_g_n_p - vs_good_pat) / sig_i ) / ( exp( (v_sellToiBuyer_g_n_p - vs_good_pat) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_g_n_p_b = exp( (v_sellToiBuyer_g_n_p - vs_bad_pat ) / sig_i ) / ( exp( (v_sellToiBuyer_g_n_p - vs_bad_pat ) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_g_n_i_g = exp( (v_sellToiBuyer_g_n_i - vs_good_imp) / sig_i ) / ( exp( (v_sellToiBuyer_g_n_i - vs_good_imp) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_g_n_i_b = exp( (v_sellToiBuyer_g_n_i - vs_bad_imp ) / sig_i ) / ( exp( (v_sellToiBuyer_g_n_i - vs_bad_imp ) / sig_i ) + 1 )    
      
      pi_sell_to_ibuyer_b_p_p_g = exp( (v_sellToiBuyer_b_p_p - vs_good_pat) / sig_i ) / ( exp( (v_sellToiBuyer_b_p_p - vs_good_pat) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_b_p_p_b = exp( (v_sellToiBuyer_b_p_p - vs_bad_pat ) / sig_i ) / ( exp( (v_sellToiBuyer_b_p_p - vs_bad_pat ) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_b_p_i_g = exp( (v_sellToiBuyer_b_p_i - vs_good_imp) / sig_i ) / ( exp( (v_sellToiBuyer_b_p_i - vs_good_imp) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_b_p_i_b = exp( (v_sellToiBuyer_b_p_i - vs_bad_imp ) / sig_i ) / ( exp( (v_sellToiBuyer_b_p_i - vs_bad_imp ) / sig_i ) + 1 )    
      
      pi_sell_to_ibuyer_b_n_p_g = exp( (v_sellToiBuyer_b_n_p - vs_good_pat) / sig_i ) / ( exp( (v_sellToiBuyer_b_n_p - vs_good_pat) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_b_n_p_b = exp( (v_sellToiBuyer_b_n_p - vs_bad_pat ) / sig_i ) / ( exp( (v_sellToiBuyer_b_n_p - vs_bad_pat ) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_b_n_i_g = exp( (v_sellToiBuyer_b_n_i - vs_good_imp) / sig_i ) / ( exp( (v_sellToiBuyer_b_n_i - vs_good_imp) / sig_i ) + 1 )    
      pi_sell_to_ibuyer_b_n_i_b = exp( (v_sellToiBuyer_b_n_i - vs_bad_imp ) / sig_i ) / ( exp( (v_sellToiBuyer_b_n_i - vs_bad_imp ) / sig_i ) + 1 )    
      
       
      
      # Laws of motion
      # 1. Calculate "inflow" rate for sellers.
      inflow_p_g = pars$mu * (  m_h_g * (1 - p.imp) * (1 - p_switch_types) * ( (1 - signal_noise) * (1 - pi_sell_to_ibuyer_g_p_p_g) + (    signal_noise) * (1 - pi_sell_to_ibuyer_g_n_p_g) ) +  # People who had good houses |--> Good & Patient
                                m_h_b * (1 - p.imp) * (    p_switch_types) * ( (1 - signal_noise) * (1 - pi_sell_to_ibuyer_b_p_p_g) + (    signal_noise) * (1 - pi_sell_to_ibuyer_b_n_p_g) ))   # People who had bad houses  |--> Good & Patient                    
      inflow_i_g = pars$mu * (  m_h_g * (    p.imp) * (1 - p_switch_types) * ( (1 - signal_noise) * (1 - pi_sell_to_ibuyer_g_p_i_g) + (    signal_noise) * (1 - pi_sell_to_ibuyer_g_n_i_g) ) +  # People who had good houses |--> Good & Impatient
                                m_h_b * (    p.imp) * (    p_switch_types) * ( (1 - signal_noise) * (1 - pi_sell_to_ibuyer_b_p_i_g) + (    signal_noise) * (1 - pi_sell_to_ibuyer_b_n_i_g) ))   # People who had bad houses  |--> Good & Impatient                         
      inflow_p_b = pars$mu * (  m_h_g * (1 - p.imp) * (    p_switch_types) * ( (    signal_noise) * (1 - pi_sell_to_ibuyer_g_p_p_b) + (1 - signal_noise) * (1 - pi_sell_to_ibuyer_g_n_p_b) ) +  # People who had good houses |--> Bad  & Patient
                                m_h_b * (1 - p.imp) * (1 - p_switch_types) * ( (    signal_noise) * (1 - pi_sell_to_ibuyer_b_p_p_b) + (1 - signal_noise) * (1 - pi_sell_to_ibuyer_b_n_p_b) ))   # People who had bad houses  |--> Bad  & Patient                    
      inflow_i_b = pars$mu * (  m_h_g * (    p.imp) * (    p_switch_types) * ( (    signal_noise) * (1 - pi_sell_to_ibuyer_g_p_i_b) + (1 - signal_noise) * (1 - pi_sell_to_ibuyer_g_n_i_b) ) +  # People who had good houses |--> Bad  & Impatient
                                m_h_b * (    p.imp) * (1 - p_switch_types) * ( (    signal_noise) * (1 - pi_sell_to_ibuyer_b_p_i_b) + (1 - signal_noise) * (1 - pi_sell_to_ibuyer_b_n_i_b) ))   # People who had bad houses  |--> Bad  & Impatient                         
      
      # 2. Calculate number of buyers directly from (buyers = sellers) and (fixed population size)
      m_b    = (pars$M - m_h_g - m_h_b + m_ib_g_e + m_ib_g_l + m_ib_b_e + m_ib_b_l)/2
      
      # 3. Calculate masses of sellers from seller flow equations + masses of buyers
      m_s_pg = inflow_p_g / (pars$lambda * m_b * pi_buy_from_good_pat)
      m_s_ig = inflow_i_g / (pars$lambda * m_b * pi_buy_from_good_imp) 
      m_s_pb = inflow_p_b / (pars$lambda * m_b * pi_buy_from_bad_pat ) 
      m_s_ib = inflow_i_b / (pars$lambda * m_b * pi_buy_from_bad_imp ) 
        
      # 4. iBuyer Inflows
      inflow_ib_g_e = pars$mu * (  m_h_g * (1 - p.imp) * (1 - p_switch_types) * ( (1 - signal_noise) * (    pi_sell_to_ibuyer_g_p_p_g) + (    signal_noise) * (    pi_sell_to_ibuyer_g_n_p_g) ) +  # People who had good houses |--> Good & Patient
                                   m_h_b * (1 - p.imp) * (    p_switch_types) * ( (1 - signal_noise) * (    pi_sell_to_ibuyer_b_p_p_g) + (    signal_noise) * (    pi_sell_to_ibuyer_b_n_p_g) ) +  # People who had bad houses  |--> Good & Patient                    
                                   m_h_g * (    p.imp) * (1 - p_switch_types) * ( (1 - signal_noise) * (    pi_sell_to_ibuyer_g_p_i_g) + (    signal_noise) * (    pi_sell_to_ibuyer_g_n_i_g) ) +  # People who had good houses |--> Good & Impatient
                                   m_h_b * (    p.imp) * (    p_switch_types) * ( (1 - signal_noise) * (    pi_sell_to_ibuyer_b_p_i_g) + (    signal_noise) * (    pi_sell_to_ibuyer_b_n_i_g) )) 
      
      inflow_ib_b_e = pars$mu * (  m_h_g * (1 - p.imp) * (    p_switch_types) * ( (    signal_noise) * (    pi_sell_to_ibuyer_g_p_p_b) + (1 - signal_noise) * (    pi_sell_to_ibuyer_g_n_p_b) ) +  # People who had good houses |--> Bad  & Patient
                                   m_h_b * (1 - p.imp) * (1 - p_switch_types) * ( (    signal_noise) * (    pi_sell_to_ibuyer_b_p_p_b) + (1 - signal_noise) * (    pi_sell_to_ibuyer_b_n_p_b) ) +  # People who had bad houses  |--> Bad  & Patient                    
                                   m_h_g * (    p.imp) * (    p_switch_types) * ( (    signal_noise) * (    pi_sell_to_ibuyer_g_p_i_b) + (1 - signal_noise) * (    pi_sell_to_ibuyer_g_n_i_b) ) +  # People who had good houses |--> Bad  & Impatient
                                   m_h_b * (    p.imp) * (1 - p_switch_types) * ( (    signal_noise) * (    pi_sell_to_ibuyer_b_p_i_b) + (1 - signal_noise) * (    pi_sell_to_ibuyer_b_n_i_b) ))   # People who had bad houses  |--> Bad  & Impatient                         
      
      inflow_ib_g_l = pars$nu * m_ib_g_e
      inflow_ib_b_l = pars$nu * m_ib_b_e
      
      
      
      # 5. iBuyer Outflows
      outflow_ib_g_e = pars$lambda * pars$lambda_ib * m_ib_g_e * m_b * pi_buy_from_good_ib_e + pars$nu * m_ib_g_e 
      outflow_ib_b_e = pars$lambda * pars$lambda_ib * m_ib_b_e * m_b * pi_buy_from_bad_ib_e  + pars$nu * m_ib_b_e 
      outflow_ib_g_l = pars$lambda * pars$lambda_ib * m_ib_g_l * m_b * pi_buy_from_good_ib_l 
      outflow_ib_b_l = pars$lambda * pars$lambda_ib * m_ib_b_l * m_b * pi_buy_from_bad_ib_l  
      
      # Calculate buyers' values from meeting the 8 seller types
      Eu_meet_hh_g_p = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_good_pat ))
      Eu_meet_hh_g_i = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_good_imp )) 
      Eu_meet_hh_b_p = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_bad_pat  )) 
      Eu_meet_hh_b_i = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_bad_imp  )) 
      Eu_meet_ib_g_e = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_good_ib_e)) 
      Eu_meet_ib_g_l = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_good_ib_l)) 
      Eu_meet_ib_b_e = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_bad_ib_e )) 
      Eu_meet_ib_b_l = sig_m * (.5772 - pars$k) + log(1 + exp(v_buy_from_bad_ib_l ))
      
      # Calculate iBuyer value functions
      
      # IBUYER_OVERHEAD = 0.02 of house price.
      
      IBUYER_OVERHEAD = pars$ibuyer.overhead
      
      
      v_ib_g_l = (pars$u_ibuyer_shock + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_good_ib_l * (p_ib_s_g_l * (1-IBUYER_OVERHEAD) )                     ) / (pars$rho + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_good_ib_l          )
      v_ib_g_e = (pars$u_ibuyer       + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_good_ib_e * (p_ib_s_g_e * (1-IBUYER_OVERHEAD) ) + pars$nu * v_ib_g_l) / (pars$rho + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_good_ib_e + pars$nu)
      v_ib_b_l = (pars$u_ibuyer_shock + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_bad_ib_l  * (p_ib_s_b_l * (1-IBUYER_OVERHEAD) )                     ) / (pars$rho + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_bad_ib_l           ) 
      v_ib_b_e = (pars$u_ibuyer       + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_bad_ib_e  * (p_ib_s_b_e * (1-IBUYER_OVERHEAD) ) + pars$nu * v_ib_b_l) / (pars$rho + pars$lambda_ib * pars$fs(b = m_b,s = 1) * pi_buy_from_bad_ib_e  + pars$nu) 
        
      ## Equations to set equal to zero ##
      #
      # Note: There are 22---
      #
      # (1) 12 prices, pinned down by FOC
      # (2) 5 value functions, pinned down by the value function equations
      # (3) 6 masses, pinned down by flow equations.
      
      # HH Price FOC
      eqn_foc_hh_gp     = - 1 / sig_m * (1 - pi_buy_from_good_pat)   * (vb + p_hh_good_pat - vs_good_pat)   + 1 # Patient   Homeowner with good house
      eqn_foc_hh_gi     = - 1 / sig_m * (1 - pi_buy_from_good_imp)   * (vb + p_hh_good_imp - vs_good_imp)   + 1 # Impatient Homeowner with good house
      eqn_foc_hh_bp     = - 1 / sig_m * (1 - pi_buy_from_bad_pat )   * (vb + p_hh_bad_pat  - vs_bad_pat )   + 1 # Patient   Homeowner with bad  house
      eqn_foc_hh_bi     = - 1 / sig_m * (1 - pi_buy_from_bad_imp )   * (vb + p_hh_bad_imp  - vs_bad_imp )   + 1 # Impatient Homeowner with bad  house
      
      # iBuyer (listing) Price FOC
      eqn_foc_ib_ge   = - 1 / sig_m * (1 - pi_buy_from_good_ib_e)  * (p_ib_s_g_e - v_ib_g_e) + 1                # iBuyer's FOC when listing a good house pre-shock
      eqn_foc_ib_gl   = - 1 / sig_m * (1 - pi_buy_from_good_ib_l)  * (p_ib_s_g_l - v_ib_g_l) + 1                # iBuyer's FOC when listing a good house post-shock
      eqn_foc_ib_be   = - 1 / sig_m * (1 - pi_buy_from_bad_ib_e )  * (p_ib_s_b_e - v_ib_b_e) + 1                # iBuyer's FOC when listing a bad  house pre-shock
      eqn_foc_ib_bl   = - 1 / sig_m * (1 - pi_buy_from_bad_ib_l )  * (p_ib_s_b_l - v_ib_b_l) + 1                # iBuyer's FOC when listing a bad  house post-shock
      
     
      
      # iBuyer purchasing price FOC.
      eqn_foc_g_p      = (1 - pars$p.imp) * (    pg_g_p) * ( ddf / sig_i * pi_sell_to_ibuyer_g_p_p_g * (1 - pi_sell_to_ibuyer_g_p_p_g) * (v_ib_g_e - p_ib_b_p_g) - pi_sell_to_ibuyer_g_p_p_g) +  # Positive signal, had good house; good type and patient
                         (    pars$p.imp) * (    pg_g_p) * ( ddf / sig_i * pi_sell_to_ibuyer_g_p_i_g * (1 - pi_sell_to_ibuyer_g_p_i_g) * (v_ib_g_e - p_ib_b_p_g) - pi_sell_to_ibuyer_g_p_i_g) +  # Positive signal, had good house; good type and impatient
                         (1 - pars$p.imp) * (1 - pg_g_p) * ( ddf / sig_i * pi_sell_to_ibuyer_g_p_p_b * (1 - pi_sell_to_ibuyer_g_p_p_b) * (v_ib_b_e - p_ib_b_p_g) - pi_sell_to_ibuyer_g_p_p_b) +  # Positive signal, had good house; bad  type and patient
                         (    pars$p.imp) * (1 - pg_g_p) * ( ddf / sig_i * pi_sell_to_ibuyer_g_p_i_b * (1 - pi_sell_to_ibuyer_g_p_i_b) * (v_ib_b_e - p_ib_b_p_g) - pi_sell_to_ibuyer_g_p_i_b)    # Positive signal, had good house; bad  type and impatient
                          
      eqn_foc_g_n      = (1 - pars$p.imp) * (    pg_g_n) * ( ddf / sig_i * pi_sell_to_ibuyer_g_n_p_g * (1 - pi_sell_to_ibuyer_g_n_p_g) * (v_ib_g_e - p_ib_b_n_g) - pi_sell_to_ibuyer_g_n_p_g) +  # Negative signal, had good house; good type and patient
                         (    pars$p.imp) * (    pg_g_n) * ( ddf / sig_i * pi_sell_to_ibuyer_g_n_i_g * (1 - pi_sell_to_ibuyer_g_n_i_g) * (v_ib_g_e - p_ib_b_n_g) - pi_sell_to_ibuyer_g_n_i_g) +  # Negative signal, had good house; good type and impatient
                         (1 - pars$p.imp) * (1 - pg_g_n) * ( ddf / sig_i * pi_sell_to_ibuyer_g_n_p_b * (1 - pi_sell_to_ibuyer_g_n_p_b) * (v_ib_b_e - p_ib_b_n_g) - pi_sell_to_ibuyer_g_n_p_b) +  # Negative signal, had good house; bad  type and patient
                         (    pars$p.imp) * (1 - pg_g_n) * ( ddf / sig_i * pi_sell_to_ibuyer_g_n_i_b * (1 - pi_sell_to_ibuyer_g_n_i_b) * (v_ib_b_e - p_ib_b_n_g) - pi_sell_to_ibuyer_g_n_i_b)    # Negative signal, had good house; bad  type and impatient
      
      eqn_foc_b_p      = (1 - pars$p.imp) * (    pg_b_p) * ( ddf / sig_i * pi_sell_to_ibuyer_b_p_p_g * (1 - pi_sell_to_ibuyer_b_p_p_g) * (v_ib_g_e - p_ib_b_p_b) - pi_sell_to_ibuyer_b_p_p_g) +  # Positive signal, had bad  house; good type and patient
                         (    pars$p.imp) * (    pg_b_p) * ( ddf / sig_i * pi_sell_to_ibuyer_b_p_i_g * (1 - pi_sell_to_ibuyer_b_p_i_g) * (v_ib_g_e - p_ib_b_p_b) - pi_sell_to_ibuyer_b_p_i_g) +  # Positive signal, had bad  house; good type and impatient
                         (1 - pars$p.imp) * (1 - pg_b_p) * ( ddf / sig_i * pi_sell_to_ibuyer_b_p_p_b * (1 - pi_sell_to_ibuyer_b_p_p_b) * (v_ib_b_e - p_ib_b_p_b) - pi_sell_to_ibuyer_b_p_p_b) +  # Positive signal, had bad  house; bad  type and patient
                         (    pars$p.imp) * (1 - pg_b_p) * ( ddf / sig_i * pi_sell_to_ibuyer_b_p_i_b * (1 - pi_sell_to_ibuyer_b_p_i_b) * (v_ib_b_e - p_ib_b_p_b) - pi_sell_to_ibuyer_b_p_i_b)    # Positive signal, had bad  house; bad  type and impatient
      
      eqn_foc_b_n      = (1 - pars$p.imp) * (    pg_b_n) * ( ddf / sig_i * pi_sell_to_ibuyer_b_n_p_g * (1 - pi_sell_to_ibuyer_b_n_p_g) * (v_ib_g_e - p_ib_b_n_b) - pi_sell_to_ibuyer_b_n_p_g) +  # Negative signal, had bad  house; good type and patient
                         (    pars$p.imp) * (    pg_b_n) * ( ddf / sig_i * pi_sell_to_ibuyer_b_n_i_g * (1 - pi_sell_to_ibuyer_b_n_i_g) * (v_ib_g_e - p_ib_b_n_b) - pi_sell_to_ibuyer_b_n_i_g) +  # Negative signal, had bad  house; good type and impatient
                         (1 - pars$p.imp) * (1 - pg_b_n) * ( ddf / sig_i * pi_sell_to_ibuyer_b_n_p_b * (1 - pi_sell_to_ibuyer_b_n_p_b) * (v_ib_b_e - p_ib_b_n_b) - pi_sell_to_ibuyer_b_n_p_b) +  # Negative signal, had bad  house; bad  type and patient
                         (    pars$p.imp) * (1 - pg_b_n) * ( ddf / sig_i * pi_sell_to_ibuyer_b_n_i_b * (1 - pi_sell_to_ibuyer_b_n_i_b) * (v_ib_b_e - p_ib_b_n_b) - pi_sell_to_ibuyer_b_n_i_b)    # Negative signal, had bad  house; bad  type and impatient
      
      
      
      
      # Sellers' (listers') value functions
      eqn_vs_good_pat       = -vs_good_pat             + (pars$u_unm_pat  +  pars$fs(b = m_b,s = 1)    * (pi_buy_from_good_pat * (p_hh_good_pat + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_from_good_pat) # You may find a guy to sell to :)
      eqn_vs_good_imp       = -vs_good_imp             + (pars$u_unm_imp  +  pars$fs(b = m_b,s = 1)    * (pi_buy_from_good_imp * (p_hh_good_imp + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_from_good_imp) # You may find a guy to sell to :)
      eqn_vs_bad_pat        = -vs_bad_pat              + (pars$u_unm_pat  +  pars$fs(b = m_b,s = 1)    * (pi_buy_from_bad_pat  * (p_hh_bad_pat  + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_from_bad_pat) # You may find a guy to sell to :)
      eqn_vs_bad_imp        = -vs_bad_imp              + (pars$u_unm_imp  +  pars$fs(b = m_b,s = 1)    * (pi_buy_from_bad_imp  * (p_hh_bad_imp  + vb)) ) / (pars$rho + pars$fs(b = m_b,s = 1) * pi_buy_from_bad_imp) # You may find a guy to sell to :)
      
      
      # Aggregate population size and housing stock balance
      eqn_buyers_equal_sellers = m_b - m_s_pg - m_s_pb - m_s_ig - m_s_ib - m_ib_g_e - m_ib_g_l - m_ib_b_e - m_ib_b_l
      eqn_good_bad_houses      = (m_h_g + m_s_pg + m_s_ig + m_ib_g_e + m_ib_g_l) - (m_h_b + m_s_pb + m_s_ib + m_ib_b_e + m_ib_b_l)

      # Laws of motion for iBuyers
      eqn_flow_ib_g_e = inflow_ib_g_e - outflow_ib_g_e
      eqn_flow_ib_b_e = inflow_ib_b_e - outflow_ib_b_e
      eqn_flow_ib_g_l = inflow_ib_g_l - outflow_ib_g_l
      eqn_flow_ib_b_l = inflow_ib_b_l - outflow_ib_b_l
      
      res = c(    100*c(eqn_foc_hh_gp,eqn_foc_hh_gi,eqn_foc_hh_bp,eqn_foc_hh_bi,eqn_foc_ib_ge,eqn_foc_ib_gl,eqn_foc_ib_be,eqn_foc_ib_bl,eqn_foc_g_p,eqn_foc_g_n,eqn_foc_b_p,eqn_foc_b_n), # Pricing FOC
                   c(eqn_vs_good_pat,eqn_vs_good_imp,eqn_vs_bad_pat,eqn_vs_bad_imp), # Value functions
                   c(eqn_buyers_equal_sellers,eqn_good_bad_houses,eqn_flow_ib_g_e,eqn_flow_ib_b_e,eqn_flow_ib_g_l,eqn_flow_ib_b_l)    )
    
      if(!return_stats) {
        return(res)
      }   else {
        
        stats = list()
        
        # Household prices
        flow_rate_hh_2_hh_g_p = pars$lambda * m_b * m_s_pg * pi_buy_from_good_pat
        flow_rate_hh_2_hh_g_i = pars$lambda * m_b * m_s_ig * pi_buy_from_good_imp
        flow_rate_hh_2_hh_b_p = pars$lambda * m_b * m_s_pb * pi_buy_from_bad_pat 
        flow_rate_hh_2_hh_b_i = pars$lambda * m_b * m_s_ib * pi_buy_from_bad_imp 
        hh_2_matched_flow_rate = flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_g_i + flow_rate_hh_2_hh_b_p + flow_rate_hh_2_hh_b_i
        
        
        # Diagnostics only
        stats$pi_buy_from_good_pat = pi_buy_from_good_pat
        stats$pi_buy_from_good_impd = pi_buy_from_good_imp
        stats$pi_buy_from_bad_pat = pi_buy_from_bad_pat
        stats$pi_buy_from_bad_imp = pi_buy_from_bad_imp
        
        stats$m_s_pg = m_s_pg
        stats$m_s_ig = m_s_ig
        stats$m_s_pb = m_s_pb
        stats$m_s_ib = m_s_ib
        
        stats$flow_rate_hh_2_hh_g_p = flow_rate_hh_2_hh_g_p
        stats$flow_rate_hh_2_hh_g_i = flow_rate_hh_2_hh_g_i
        stats$flow_rate_hh_2_hh_b_p = flow_rate_hh_2_hh_b_p
        stats$flow_rate_hh_2_hh_b_i = flow_rate_hh_2_hh_b_i 
        stats$hh_2_matched_flow_rate = hh_2_matched_flow_rate
          
        
        stats$price_hh_good_patient   = p_hh_good_pat 
        stats$price_hh_good_impatient = p_hh_good_imp
        stats$price_hh_bad_patient    = p_hh_bad_pat
        stats$price_hh_bad_impatient  = p_hh_bad_imp
        stats$price_hh_mean           = 1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_p * p_hh_good_pat + flow_rate_hh_2_hh_g_i * p_hh_good_imp + flow_rate_hh_2_hh_b_p * p_hh_bad_pat + flow_rate_hh_2_hh_b_i * p_hh_bad_imp)
        stats$price_hh_patient_mean   = 1 / (flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_b_p) * (flow_rate_hh_2_hh_g_p * p_hh_good_pat   + flow_rate_hh_2_hh_b_p * p_hh_bad_pat)
        stats$price_hh_impatient_mean = 1 / (flow_rate_hh_2_hh_g_i + flow_rate_hh_2_hh_b_i) * (flow_rate_hh_2_hh_g_i * p_hh_good_imp   + flow_rate_hh_2_hh_b_i * p_hh_bad_imp)
        stats$price_hh_good_mean      = 1 / (flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_g_i) * (flow_rate_hh_2_hh_g_p * p_hh_good_pat   + flow_rate_hh_2_hh_g_i * p_hh_good_imp)
        stats$price_hh_bad_mean       = 1 / (flow_rate_hh_2_hh_b_p + flow_rate_hh_2_hh_b_i) * (flow_rate_hh_2_hh_b_p * p_hh_bad_pat    + flow_rate_hh_2_hh_b_i * p_hh_bad_imp)
        
        
        
        # iBuyer purchase prices
        flow_rate_hh_2_ib_p_g = pars$mu * m_h_g * ( (1 - p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_g + (    p.imp) * pi_sell_to_ibuyer_g_p_i_g)  + 
                                                    (    p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_b + (    p.imp) * pi_sell_to_ibuyer_g_p_i_b))
        flow_rate_hh_2_ib_p_b = pars$mu * m_h_b * ( (1 - p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_b + (    p.imp) * pi_sell_to_ibuyer_b_p_i_b)  + 
                                                    (    p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_g + (    p.imp) * pi_sell_to_ibuyer_b_p_i_g))
        flow_rate_hh_2_ib_n_g = pars$mu * m_h_g * ( (1 - p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_g + (    p.imp) * pi_sell_to_ibuyer_g_n_i_g)  + 
                                                    (    p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_b + (    p.imp) * pi_sell_to_ibuyer_g_n_i_b))
        flow_rate_hh_2_ib_n_b = pars$mu * m_h_b * ( (1 - p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_b + (    p.imp) * pi_sell_to_ibuyer_b_n_i_b)  + 
                                                    (    p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_g + (    p.imp) * pi_sell_to_ibuyer_b_n_i_g))
        
        seller_inflow = (m_h_g + m_h_b) * pars$mu
        ibuyer_inflow = inflow_ib_g_e + inflow_ib_b_e
        pi_sell_to_ibuyer = ibuyer_inflow / seller_inflow 
        stats$ib_share = pi_sell_to_ibuyer
        
        stats$price_ib_buy_p_g  = p_ib_b_p_g
        stats$price_ib_buy_p_b  = p_ib_b_p_b
        stats$price_ib_buy_n_g  = p_ib_b_n_g
        stats$price_ib_buy_n_b  = p_ib_b_n_b
        stats$price_ib_buy_mean = 1 / ibuyer_inflow * (flow_rate_hh_2_ib_p_g * p_ib_b_p_g + flow_rate_hh_2_ib_p_b * p_ib_b_p_b + flow_rate_hh_2_ib_n_g * p_ib_b_n_g + flow_rate_hh_2_ib_n_b * p_ib_b_n_b)

        
        
        # Calculate average prices that iBuyers are paying for good and bad houses. 
        flow_rate_good_house_sold_to_ibuyer = pars$mu * m_h_g * ( (1 - p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_g + (    p.imp) * pi_sell_to_ibuyer_g_p_i_g)) +
                                              pars$mu * m_h_b * ( (    p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_g + (    p.imp) * pi_sell_to_ibuyer_b_p_i_g)) + 
                                              pars$mu * m_h_g * ( (1 - p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_g + (    p.imp) * pi_sell_to_ibuyer_g_n_i_g)) + 
                                              pars$mu * m_h_b * ( (    p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_g + (    p.imp) * pi_sell_to_ibuyer_b_n_i_g))
          
        flow_rate_bad_house_sold_to_ibuyer = pars$mu * m_h_g * ( (    p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_b + (    p.imp) * pi_sell_to_ibuyer_g_p_i_b)) +
                                             pars$mu * m_h_b * ( (1 - p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_b + (    p.imp) * pi_sell_to_ibuyer_b_p_i_b)) + 
                                             pars$mu * m_h_g * ( (    p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_b + (    p.imp) * pi_sell_to_ibuyer_g_n_i_b)) + 
                                             pars$mu * m_h_b * ( (1 - p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_b + (    p.imp) * pi_sell_to_ibuyer_b_n_i_b))

        
        stats$price_ib_buy_good_mean =    (pars$mu * m_h_g * ( (1 - p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_g + (    p.imp) * pi_sell_to_ibuyer_g_p_i_g)) * p_ib_b_p_g +
                                           pars$mu * m_h_b * ( (    p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_g + (    p.imp) * pi_sell_to_ibuyer_b_p_i_g)) * p_ib_b_p_b + 
                                           pars$mu * m_h_g * ( (1 - p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_g + (    p.imp) * pi_sell_to_ibuyer_g_n_i_g)) * p_ib_b_n_g + 
                                           pars$mu * m_h_b * ( (    p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_g + (    p.imp) * pi_sell_to_ibuyer_b_n_i_g)) * p_ib_b_n_b) / flow_rate_good_house_sold_to_ibuyer
        stats$price_ib_buy_bad_mean  =    (pars$mu * m_h_g * ( (    p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_b + (    p.imp) * pi_sell_to_ibuyer_g_p_i_b)) * p_ib_b_p_g +
                                           pars$mu * m_h_b * ( (1 - p_switch_types) * (    signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_b + (    p.imp) * pi_sell_to_ibuyer_b_p_i_b)) * p_ib_b_p_b + 
                                           pars$mu * m_h_g * ( (    p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_b + (    p.imp) * pi_sell_to_ibuyer_g_n_i_b)) * p_ib_b_n_g + 
                                           pars$mu * m_h_b * ( (1 - p_switch_types) * (1 - signal_noise) * ( (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_b + (    p.imp) * pi_sell_to_ibuyer_b_n_i_b)) * p_ib_b_n_b)  / flow_rate_bad_house_sold_to_ibuyer
                                        
        
        # To calculate iBuyer buy discount, we  want:
        # (price ibuyer pays | observables) / (hh price | observables)
        #
        # Price iBuyer pays | observables is easy (it's direct)
        #
        # (hh price | observables) is harder
        hh_price_given_g_p = (p_hh_good_pat * (1 - p_switch_types) * (1 - signal_noise) * (1 - p.imp) +
                              p_hh_good_imp * (1 - p_switch_types) * (1 - signal_noise) * (    p.imp) + 
                              p_hh_bad_pat  * (    p_switch_types) * (    signal_noise) * (1 - p.imp) +
                              p_hh_bad_imp  * (    p_switch_types) * (    signal_noise) * (    p.imp)  ) / ( (1 - p_switch_types) * (1 - signal_noise) + (    p_switch_types) * (    signal_noise))
        
        hh_price_given_g_n = (p_hh_good_pat * (1 - p_switch_types) * (    signal_noise) * (1 - p.imp) +
                              p_hh_good_imp * (1 - p_switch_types) * (    signal_noise) * (    p.imp) + 
                              p_hh_bad_pat  * (    p_switch_types) * (1 - signal_noise) * (1 - p.imp) +
                              p_hh_bad_imp  * (    p_switch_types) * (1 - signal_noise) * (    p.imp)  ) / ( (1 - p_switch_types) * (    signal_noise) + (    p_switch_types) * (1 - signal_noise))
        
        hh_price_given_b_p = (p_hh_good_pat * (    p_switch_types) * (1 - signal_noise) * (1 - p.imp) +
                              p_hh_good_imp * (    p_switch_types) * (1 - signal_noise) * (    p.imp) + 
                              p_hh_bad_pat  * (1 - p_switch_types) * (    signal_noise) * (1 - p.imp) +
                              p_hh_bad_imp  * (1 - p_switch_types) * (    signal_noise) * (    p.imp)  ) / ( (    p_switch_types) * (1 - signal_noise) + (1 - p_switch_types) * (    signal_noise))
        
        hh_price_given_b_n = (p_hh_good_pat * (    p_switch_types) * (    signal_noise) * (1 - p.imp) +
                              p_hh_good_imp * (    p_switch_types) * (    signal_noise) * (    p.imp) + 
                              p_hh_bad_pat  * (1 - p_switch_types) * (1 - signal_noise) * (1 - p.imp) +
                              p_hh_bad_imp  * (1 - p_switch_types) * (1 - signal_noise) * (    p.imp)  ) / ( (    p_switch_types) * (    signal_noise) + (1 - p_switch_types) * (1 - signal_noise))
        
          
        ib_discount_gp = p_ib_b_p_g / hh_price_given_g_p - 1
        ib_discount_gn = p_ib_b_n_g / hh_price_given_g_n - 1
        ib_discount_bp = p_ib_b_p_b / hh_price_given_b_p - 1
        ib_discount_bn = p_ib_b_n_b / hh_price_given_b_n - 1
                                        
        
        
        
        # iBuyer sale prices
        flow_rate_ib_2_hh_g_e = pars$lambda * pars$lambda_ib * m_ib_g_e * m_b * pi_buy_from_good_ib_e
        flow_rate_ib_2_hh_b_e = pars$lambda * pars$lambda_ib * m_ib_b_e * m_b * pi_buy_from_bad_ib_e
        flow_rate_ib_2_hh_g_l = pars$lambda * pars$lambda_ib * m_ib_g_l * m_b * pi_buy_from_good_ib_l
        flow_rate_ib_2_hh_b_l = pars$lambda * pars$lambda_ib * m_ib_b_l * m_b * pi_buy_from_bad_ib_l
        ibuyer_outflow = flow_rate_ib_2_hh_g_e + flow_rate_ib_2_hh_b_e + flow_rate_ib_2_hh_g_l + flow_rate_ib_2_hh_b_l
        
        
        stats$price_ib_sell_g_e = p_ib_s_g_e
        stats$price_ib_sell_b_e = p_ib_s_b_e
        stats$price_ib_sell_g_l = p_ib_s_g_l
        stats$price_ib_sell_b_l = p_ib_s_b_l
        stats$price_ib_sell_good_mean = 1 / (flow_rate_ib_2_hh_g_e + flow_rate_ib_2_hh_g_l) * (flow_rate_ib_2_hh_g_e * p_ib_s_g_e + flow_rate_ib_2_hh_g_l * p_ib_s_g_l)
        stats$price_ib_sell_bad_mean = 1 /  (flow_rate_ib_2_hh_b_e + flow_rate_ib_2_hh_b_l) * (flow_rate_ib_2_hh_b_e * p_ib_s_b_e + flow_rate_ib_2_hh_b_l * p_ib_s_b_l)
        stats$price_ib_sell_mean = 1 / ibuyer_outflow * (flow_rate_ib_2_hh_g_e * p_ib_s_g_e + flow_rate_ib_2_hh_b_e * p_ib_s_b_e + flow_rate_ib_2_hh_g_l * p_ib_s_g_l + flow_rate_ib_2_hh_b_l * p_ib_s_b_l)
    
        ibuyer_sell_premium_good = stats$price_ib_sell_good_mean / stats$price_hh_good_mean - 1
        ibuyer_sell_premium_bad  = stats$price_ib_sell_bad_mean  / stats$price_hh_bad_mean  - 1
        
        
        ## AVERAGE CONDITIONAL PREMIUM CHARGED
        stats$ibuyer_sell_premium_AVERAGE = (ibuyer_sell_premium_good * flow_rate_good_house_sold_to_ibuyer + ibuyer_sell_premium_bad * flow_rate_bad_house_sold_to_ibuyer) / (flow_rate_good_house_sold_to_ibuyer + flow_rate_bad_house_sold_to_ibuyer)
        
        ## Value functions
        stats$vs_good_pat = vs_good_pat
        stats$vs_good_imp = vs_good_imp
        stats$vs_bad_pat  = vs_bad_pat
        stats$vs_bad_imp  = vs_bad_imp
        stats$vb          = vb
        stats$vh_good     = v_h_g
        stats$vh_bad      = v_h_b
        stats$vhh         = m_h_g * v_h_g + m_h_b * v_h_b + m_b * vb + m_s_pg * vs_good_pat + m_s_pb * vs_bad_pat + m_s_ig * vs_good_imp + m_s_ib * vs_bad_imp
        
        
        
        # HH Time on Market
        stats$tom_hh_good_patient     = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_from_good_pat) * 365
        stats$tom_hh_good_impatient   = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_from_good_imp) * 365 
        stats$tom_hh_bad_patient      = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_from_bad_pat ) * 365 
        stats$tom_hh_bad_impatient    = 1 / (pars$fs(b = m_b,s = 1)  * pi_buy_from_bad_imp ) * 365 
        stats$tom_hh_mean             = 1 / hh_2_matched_flow_rate * (flow_rate_hh_2_hh_g_i * stats$tom_hh_good_impatient + flow_rate_hh_2_hh_b_i * stats$tom_hh_bad_impatient + flow_rate_hh_2_hh_g_p * stats$tom_hh_good_patient + flow_rate_hh_2_hh_b_p * stats$tom_hh_bad_patient) 
        stats$tom_hh_patient_mean     = 1 / (flow_rate_hh_2_hh_g_p + flow_rate_hh_2_hh_b_p) * (flow_rate_hh_2_hh_g_p * stats$tom_hh_good_patient   + flow_rate_hh_2_hh_b_p * stats$tom_hh_bad_patient  )
        stats$tom_hh_impatient_mean   = 1 / (flow_rate_hh_2_hh_g_i + flow_rate_hh_2_hh_b_i) * (flow_rate_hh_2_hh_g_i * stats$tom_hh_good_impatient + flow_rate_hh_2_hh_b_i * stats$tom_hh_bad_impatient)
        
        flow_rate_hh_2_hh_g_p = pars$lambda * m_b * m_s_pg * pi_buy_from_good_pat
        flow_rate_hh_2_hh_g_i = pars$lambda * m_b * m_s_ig * pi_buy_from_good_imp
        flow_rate_hh_2_hh_b_p = pars$lambda * m_b * m_s_pb * pi_buy_from_bad_pat 
        flow_rate_hh_2_hh_b_i = pars$lambda * m_b * m_s_ib * pi_buy_from_bad_imp 
        
        # IB Time on Market
        # TOM | no switch * p(no switch) + TOM | switch * P(switch)
        
        stats$tom_ib_g_e   = 1 / (pars$lambda_ib * pars$fs(b = m_b,s = 1)  * pi_buy_from_good_ib_e) * 365
        stats$tom_ib_g_l   = 1 / (pars$lambda_ib * pars$fs(b = m_b,s = 1)  * pi_buy_from_good_ib_l) * 365 
        stats$tom_ib_b_e   = 1 / (pars$lambda_ib * pars$fs(b = m_b,s = 1)  * pi_buy_from_bad_ib_e ) * 365 
        stats$tom_ib_b_l   = 1 / (pars$lambda_ib * pars$fs(b = m_b,s = 1)  * pi_buy_from_bad_ib_l ) * 365 
        
        
        stats$v_ib_g_e = v_ib_g_e
        stats$v_ib_b_e = v_ib_b_e
        stats$v_ib_mean = 1 / (v_ib_g_e + v_ib_b_e) * (stats$v_ib_g_e * flow_rate_ib_2_hh_g_e + stats$v_ib_b_e * flow_rate_ib_2_hh_b_e)
          
        
        
        stats$ib_flow_quantity       = ibuyer_outflow
        stats$ib_flow_quantity_early = flow_rate_ib_2_hh_g_e + flow_rate_ib_2_hh_b_e
        stats$ib_flow_quantity_late  = flow_rate_ib_2_hh_g_l + flow_rate_ib_2_hh_b_l
        
        # P that IB drops prices---because it's steady state it's just the fraction of ib sales coming from late sales
        stats$ib_p_reduce_prices =  (flow_rate_ib_2_hh_g_l + flow_rate_ib_2_hh_b_l) / ibuyer_outflow
        
        
        # Total time on market is: E|no switch * P(no switch) + E|switch * P(switch)
        # E|no switch = stats$tom_ib_g_e
        # P(no switch) = (1 - stats$ib_p_reduce_prices)
        # E|switch = stats$tom_ib_g_l + 1 / pars$nu
        stats$tom_ib_g_mean = stats$tom_ib_g_e * (1 - stats$ib_p_reduce_prices) + (1 / pars$nu + stats$tom_ib_g_l) * stats$ib_p_reduce_prices
        stats$tom_ib_b_mean = stats$tom_ib_b_e * (1 - stats$ib_p_reduce_prices) + (1 / pars$nu + stats$tom_ib_b_l) * stats$ib_p_reduce_prices 
        
        stats$tom_ib_mean             = 1 / (flow_rate_ib_2_hh_g_e + flow_rate_ib_2_hh_b_e) * (stats$tom_ib_g_mean * flow_rate_ib_2_hh_g_e + stats$tom_ib_b_mean * flow_rate_ib_2_hh_b_e)
        
        # Calculate some ibuyer stats:
        stats$ibuyer.share 
        
        
        
        stats$ib_prod_share = rbind(
                                    data.table(past = 'g',signal = 'p',patient = 'p',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_g * (1 - p_switch_types) * (1 - signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_p_g), 
                                    data.table(past = 'g',signal = 'p',patient = 'i',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_g * (1 - p_switch_types) * (1 - signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_g_p_i_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_p_g),
                                    data.table(past = 'g',signal = 'p',patient = 'p',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_g * (    p_switch_types) * (    signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_g_p_p_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_p_g),
                                    data.table(past = 'g',signal = 'p',patient = 'i',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_g * (    p_switch_types) * (    signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_g_p_i_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_p_g),
        
                                    data.table(past = 'g',signal = 'n',patient = 'p',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_g * (1 - p_switch_types) * (    signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_n_g),
                                    data.table(past = 'g',signal = 'n',patient = 'i',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_g * (1 - p_switch_types) * (    signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_g_n_i_g,profitability = stats$price_ib_sell_good_mean  * (1-IBUYER_OVERHEAD)- stats$price_ib_buy_n_g),
                                    data.table(past = 'g',signal = 'n',patient = 'p',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_g * (    p_switch_types) * (1 - signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_g_n_p_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_n_g),
                                    data.table(past = 'g',signal = 'n',patient = 'i',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_g * (    p_switch_types) * (1 - signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_g_n_i_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_n_g),
          
                                    data.table(past = 'b',signal = 'p',patient = 'p',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_b * (    p_switch_types) * (1 - signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_p_b),
                                    data.table(past = 'b',signal = 'p',patient = 'i',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_b * (    p_switch_types) * (1 - signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_b_p_i_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_p_b),
                                    data.table(past = 'b',signal = 'p',patient = 'p',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_b * (1 - p_switch_types) * (    signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_b_p_p_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_p_b),
                                    data.table(past = 'b',signal = 'p',patient = 'i',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_b * (1 - p_switch_types) * (    signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_b_p_i_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_p_b),
        
                                    data.table(past = 'b',signal = 'n',patient = 'p',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_b * (    p_switch_types) * (    signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_n_b),
                                    data.table(past = 'b',signal = 'n',patient = 'i',type = 'g',share = 1/ibuyer_inflow * pars$mu * m_h_b * (    p_switch_types) * (    signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_b_n_i_g,profitability = stats$price_ib_sell_good_mean * (1-IBUYER_OVERHEAD) - stats$price_ib_buy_n_b),
                                    data.table(past = 'b',signal = 'n',patient = 'p',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_b * (1 - p_switch_types) * (1 - signal_noise) *  (1 - p.imp) * pi_sell_to_ibuyer_b_n_p_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_n_b),
                                    data.table(past = 'b',signal = 'n',patient = 'i',type = 'b',share = 1/ibuyer_inflow * pars$mu * m_h_b * (1 - p_switch_types) * (1 - signal_noise) *  (    p.imp) * pi_sell_to_ibuyer_b_n_i_b,profitability = stats$price_ib_sell_bad_mean * (1-IBUYER_OVERHEAD)  - stats$price_ib_buy_n_b)
                                )
        
        
        stats$ib_prod_share[,share_total := share * stats$ib_share]
         
          
        ## AVERAGE CONDITIONAL DISCOUNT PAID 
        stats$ibuyer_buy_discount_AVERAGE = ib_discount_gp * sum(stats$ib_prod_share[past == 'g' & signal == 'p']$share) + 
                                            ib_discount_gn * sum(stats$ib_prod_share[past == 'g' & signal == 'n']$share) + 
                                            ib_discount_bp * sum(stats$ib_prod_share[past == 'b' & signal == 'p']$share) + 
                                            ib_discount_bn * sum(stats$ib_prod_share[past == 'b' & signal == 'n']$share)
        
        
        return(stats)
        
        
      }
      
    }
    
    # Solve it.
    tol = 1e-15
    solns = list()
    
    
    # 1. If we were supplied a guess, try it.
    if(!is.null(ibuyer.guess)) {
      x0 = ibuyer.guess
      soln = nleqslv(x = x0,fn = toSolve,method = 'Broyden',global = 'gline',xscalm = 'auto',control = list(maxit = 500,xtol = 1e-25,ftol = 1e-15,stepmax = -1,allowSingular = T))
      distance = sum(soln$fvec^2)
      if(distance < tol) {
        solns = list(dist = distance,soln = soln$x)
      }
    }
    
    # Okay, that didn't work. Try a guess motivated by the no-ibuyer case
    if(length(solns) == 0) {
      x0 = c(no.ibuyer.solution[1]-0,no.ibuyer.solution[2]-0,no.ibuyer.solution[3]-0,no.ibuyer.solution[4]-0, # HH selling prices
                       no.ibuyer.solution[1]-2,no.ibuyer.solution[2]-4,no.ibuyer.solution[3]-4,no.ibuyer.solution[4]-4, # iBuyer buying prices
                       no.ibuyer.solution[1]-0,no.ibuyer.solution[2]-0,no.ibuyer.solution[3]-0,no.ibuyer.solution[4]-0, # iBuyer selling prices
                       no.ibuyer.solution[5]+1,no.ibuyer.solution[6]+1,no.ibuyer.solution[7]+1,no.ibuyer.solution[8]+1, # Value functions
                       no.ibuyer.solution[10]-.03,no.ibuyer.solution[11]-0.03,.015,.015,.01,.01) # Masses
      soln = nleqslv(x = x0,fn = toSolve,method = 'Broyden',global = 'gline',xscalm = 'auto',control = list(maxit = 250,xtol = 1e-20,ftol = 1e-10,stepmax = -1,allowSingular = T))
      distance = sum(soln$fvec^2)
      if(distance < tol) {
        solns = list(dist = distance,soln = soln$x)
      }
    }    
    
    
    # Okay, that didn't work. Yry a bunch of probabilistic iterations and pick the best one: (Very rare to get here)
    if(length(solns) == 0) { 
      solns <- foreach(ii = 1:20, .combine = c, .packages = c('nleqslv','data.table')) %dopar% {
        current.rand = rnorm(length(x0),mean=1,sd = .001) 
        this.guess   = x0 * current.rand
        soln = nleqslv(x = this.guess,fn = toSolve,method = 'Broyden',global = 'gline',xscalm = 'auto',control = list(maxit = 250,xtol = 1e-20,ftol = 1e-10,stepmax = -1,allowSingular = T))
        distance = sum(soln$fvec^2)
        if(distance < 1e-10) {
          toret = list(dist = distance,soln = soln$x)
          return(toret)
        }
      }
    }
    
    # 2. If that doesn't work, try a different algorithm
    if(length(solns) == 0) {
      solns <- foreach(ii = 1:20, .combine = c, .packages = c('nleqslv','data.table')) %dopar% {
        current.rand = rnorm(length(x0),mean=1,sd = .01) 
        this.guess   = x0 * current.rand
        soln = nleqslv(x = this.guess,fn = toSolve,method = 'Broyden',global = 'gline',xscalm = 'auto',control = list(maxit = 250,xtol = 1e-20,ftol = 1e-10,stepmax = -1,allowSingular = T))
        distance = sum(soln$fvec^2)
        if(distance < 1e-10) {
          toret = list(dist = distance,soln = soln$x)
          return(toret)
        }
      }
    }
    
    # 3. If that doesn't work, try a different algorithm
    if(length(solns) == 0) {
      solns <- foreach(ii = 1:20, .combine = c, .packages = c('nleqslv','data.table')) %dopar% {
        current.rand = rnorm(length(x0),mean=1,sd = .001) 
        this.guess   = x0 * current.rand
        soln = nleqslv(x = this.guess,fn = toSolve,method = 'Newton',global = 'gline',xscalm = 'auto',control = list(maxit = 250,xtol = 1e-20,ftol = 1e-10,stepmax = -1,allowSingular = T))
        distance = sum(soln$fvec^2)
        if(distance < 1e-10) {
          toret = list(dist = distance,soln = soln$x)
          return(toret)
        }
      }
    }
    
    # 4. If still nothing works,
    if(length(solns) == 0) {
      print('Failed to find solution to iBuyer equilibrium')
      soln = nleqslv(x = x0,fn = toSolve,method = 'Broyden',global = 'gline',xscalm = 'auto',control = list(maxit = 250,xtol = 1e-20,ftol = 1e-10,stepmax = -1,allowSingular = T))
      distance = sum(soln$fvec^2)
      solns = list(dist = distance,soln = soln$x)
    }
    
    
    # 5. Take the solved model.
    ibuyer.stats <- toSolve(solns$soln,return_stats = T)
    ibuyer.stats$DIST = solns$dist
    ibuyer.solution = solns$soln
    
    toReturn = list()
    
    toReturn$no.ibuyer$solution = no.ibuyer.solution
    toReturn$no.ibuyer$stats    = no.ibuyer.stats
    toReturn$ibuyer$solution <- ibuyer.solution
    toReturn$ibuyer$stats    <- ibuyer.stats
    
    
    return(toReturn)
  }  else {
    
    # no-ibuyer solution
    toReturn = list()
    toReturn$no.ibuyer$solution = no.ibuyer.solution
    toReturn$no.ibuyer$stats    = no.ibuyer.stats
    
    return(toReturn)
  }
  
}

calculateNoiBuyerMoments <- function(solution) {
  
  stats <- solution$no.ibuyer$stats
  
  # Calculate no ibuyer moments
  MOMENT_house.price            = stats$price_hh_mean
  MOMENT_hh.tom                 = stats$tom_hh_mean
  MOMENT_hh.tom.imp.delta       = stats$tom_hh_patient_mean - stats$tom_hh_impatient_mean
  MOMENT_hh.price.imp.delta.pct = stats$price_hh_mean_patient / stats$price_hh_mean_impatient - 1
  MOMENT_hh.price.sd.xs         = stats$log_price_hh_std_xs
  MOMENT_hh.price.sd.ts         = stats$log_price_hh_std_ts
  
  moments = rbind(                   data.table(name = 'tom.hh',              value = MOMENT_hh.tom),
                                     data.table(name = 'imp.price.delta.pct', value = MOMENT_hh.price.imp.delta.pct),
                                     data.table(name = 'imp.list.delta.days', value = MOMENT_hh.tom.imp.delta),
                                     data.table(name = 'price_listed'       , value = MOMENT_house.price),
                                     data.table(name = 'price_sd_xs'      ,   value = MOMENT_hh.price.sd.xs),
                                     data.table(name = 'price_sd_ts'      ,   value = MOMENT_hh.price.sd.ts)) 
  return(moments)
  
}

calculateiBuyerMoments <- function(solution,derivatives) {
  
  pre.ibuyer.moments <- calculateNoiBuyerMoments(solution)
  
  
  stats    <- solution$ibuyer$stats
  MOMENT_ib.share        <- stats$ib_share
  MOMENT_ib.buy.discount <- stats$ibuyer_buy_discount_AVERAGE # stats$price_ib_buy_mean  / stats$price_hh_mean - 1
  MOMENT_ib.sell.premium <- stats$ibuyer_sell_premium_AVERAGE #stats$price_ib_sell_mean / stats$price_hh_mean - 1
  MOMENT_ib.tom.delta    <- stats$tom_ib_mean - stats$tom_hh_mean
  MOMENT_ib.p.drop       <- stats$ib_p_reduce_prices
  MOMENT_ib.dp_dnoise    <- derivatives$ds_dn
  MOMENT_ib.dp_dspeed    <- derivatives$ds_dv
  MOMENT_dist            <- stats$DIST
  MOMENT_share_unprofit  <- sum(stats$ib_prod_share[profitability < 0]$share)
  
  profitability   = sum(stats$ib_prod_share$share * stats$ib_prod_share$profitability)/stats$price_hh_mean
  
  
  ibuyer.moments = rbind(            data.table(name = 'p.reduce.prices',      value = MOMENT_ib.p.drop),
                                     data.table(name = 'tom.ib.delta',         value = MOMENT_ib.tom.delta),
                                     data.table(name = 'iBuyer.share',         value = MOMENT_ib.share),
                                     data.table(name = 'buy.discount'       ,  value = MOMENT_ib.buy.discount),
                                     data.table(name = 'sell.premium'      ,   value = MOMENT_ib.sell.premium),
                                     data.table(name = 'diB.dLiq'      ,       value = MOMENT_ib.dp_dspeed),
                                     data.table(name = 'diB.dErr'      ,       value = MOMENT_ib.dp_dnoise),
                                     data.table(name = 'DIST'          ,       value = MOMENT_dist),
                                     data.table(name = 'derr.sign'     ,       value = sign(MOMENT_ib.dp_dnoise)),
                                     data.table(name = 'dliq.sign'     ,       value = sign(MOMENT_ib.dp_dspeed)),
                                     data.table(name = 'share.unprofitable',   value = MOMENT_share_unprofit),
                                     data.table(name = 'ibuyer.spread'     ,   value = profitability)) 
  
  

  return(rbind(ibuyer.moments,pre.ibuyer.moments))
  

  
}


run_all()



